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I; INTRODUCTION - • 

In recent years tremendous advances in the practice of astronomy have become 
possible by escaping the atmospheric shield surrounding earth-based astronomical 
observatories. Earth-orbiting telescopes are unaffected by the atmosphere and 
thus capable of achieving resolving power that is completely beyond the reach of 
terrestrial instruments. Refraction anomalies and background glow, heretofore 
inherent limiting factors in celestial observations, are problem areas that have 
been successfully eliminated. They have been replaced by problems associated 
with maintaining extremely low optical tolerances in equipment exposed to the 
harsh envi ronment of rocket launch and orbit. Nonetheless, the problem of main- 
taining small tolerances in construction and alignment is a solvable problem, 
whereas the problems of accuracy in land based telescopes can only be solved up 
to the limit imposed by irregular atmospheric refraction. 

As stated, the most critical problem of orbiting telescopes is the main^ 
tenance of accuracy in construction and alignment. This problem is compounded 
by the uncertainty of the disturbances which affect the telescope during launch 
and in orbit. For this reason the critical parts of the instruments must be 
designed to be as insensitive as possible to external disturbances. There are a 
number of ways this "passive control" of optical surfaces may be provided. 

In addition to the passive shielding method, it is necessary to consider 
slightly more elaborate techniques to actively control some of the most critical 
optical surfaces. The object of these "active" methods is not to replace the 
passive shielding, but to nullify the effects of those disturbances that cannot 
be predicted and adequately compensated by passive means. 

The most critical surface on an orbiting telescope is the primary mirror. 
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The limits on the resolving power of the instruments are imposed'by the size and 
precision of the surface.' For this reason, a' number of studies' have considered 
the "active" suppression of surface irregularities (1 , 2). Such an approach to 
the problem has proven to be entirely feasible and these studies are continuing. 
The active mirror studies which have been performed thus far - have considered the 
use of either external forces on the back of the mi rror^or segmented mirrors to 
control the shape of the mirror surface (3). It has been shown that an overall 
improvement in the surface figure may be achieved in this manner and that such 
a figure compensation can be performed automatically in a servo control mode (4). 

Having determined that active control of the mirror surface is feasible and 
desirable, attention is directed to the determination' of the best way to mechanize 
the controlled flexing of the surface. To this end such things as cost, weight, 
reliability and range of operation must be considered in the selection of the 
best means of implementation. 

The results of a study of the feasibility of an active method of surface 
error control using thermal elements are presented in this report. It is shown 
that the control effort of the thermal elements is sufficient for the purpose, 
and that such benefits as low cost, low weight and high reliability may be 
achieved in conjunction with a significant reduction in the mirror surface error 
figure . 
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II. ANALYTICAL APPROACH 


Introduction 

To accomplish the objectives of this research it was necessary to. formulate 
a thermoelastic response model of the primary mirror which' could be used to 
simulate the ; response of the structure to disturbances and to controlled thermal 
inputs. Iii addition the control strategies for the active system were developed 
through the use of the response -model . 

The features that were included in the response model were as follows: 

(a) steady state thermal response 

(b) transient thermal response 

(c) structural response corresponding to a and b 

(d) generation of the influence matrix (thermal input - surface de- 
flection output) for use with the control program. 

The features that were included in the control program were: 

(a) arbitrary control and output node selection 

( b ) free thermal expansion removal 

(c) optimal mean-squared thermal control computation 
■(d) surface error computation. 

Concepts of Discretization 

Due to the axi symmetric geometry of -the unperturbed mirror it was decided, 
after preliminary investigation, to use modal decomposition in one dimension 
(angular) and to use a two dimensional finite element discretization in the 
radial and longitudinal coordinates. Such a system had been shown to give 
highly reliable results in work done on solid rocket grains by Wilson (5) and 
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was familiar to one of the authors. It was' felt that such a method would be superior 
in many respects to a- complete three dimension finite element-discretization if the 
number of angular modes arising from unsymmetrieal disturbance and thermal control 
excitation were small. 

The angular span of the thermal patches were on the order of 10°, which results 

in a significant number of harmonics necessary to characterize their input effect. 

(Early in the work we were prepared' to use 100 modes in the, analysis) . However, 

the inverse of the thermal coefficient matrix' C which is obtained after the 

• -2 

discretization is completed contains elements that are proportional to n , where 
n is the order of the mode. Therefore, the temperature of the mirror is affected 
by any mode of a patch in proportion to the inverse square of the order of the 
mode. This has the effect of reducing the total number of modes used for each 
patch, since the higher order harmonics are attenuated to such an extent that 
their presence is inconsequential . It has been shown that a truncation after 
ten modes for a 10° patch angle results in a temperature error less than .1% at 
any point in the mirror. 

Due to the ability of this program to function with a very small number of 
radial modes, it is estimated that a savings of 90%' in total number of computations 
and 99% in storage locations is made over a full three dimensional finite element 
model . 

One disadvantage of the discretization scheme is the necessity for a large 
storage to compute a true transient response in time. For example, to compute 
the deformation of the structure at time t 1} it is necessary to compute the 
nodal temperatures of the two dimensional' finite element models corresponding to 
all angular modes , sum the temperatures , and compute the deformation. If the 
process is to be repeated at a' later time t 2 , then all of the nodal temperatures 
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for every mode must be computed and stored at time t x . 

This amount of storage was not available to us, so we programmed the transient 
loop to compute the structural response at times ti, t 2 , t'j, ... etc., starting 
each computation at t = 0. This replaces storaged needed by increased run time. 


Thermal Response Formulation 

The first stage of a thermoelastic response is the determination of the 
temperature field, T(r,e,z,t) that exists in the mirror. The temperature is 
governed by the partial differential equation 

v2T-fijpT t (i.i-i) 


with boundary conditions 

q * -kT on the bottom 
^ n 

q = -kT = 0 on the sides 
^ n 

and 


q = -kT = -e(T 4 - T 4 ) on the front 
^ n o 


for the geometry shown in Figure 1 1 -1 . 
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Figure I I -1 
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( 1 1 -3a) 


(IX -3b) 


The finite element discretization is obtained initially by designating nodes in 

the r, z plane and identifying triangular areas between node sets. The equations 

II -3 hold over any area in the r-z plane and in particular over the triangular 

finite element areas . Within each triangular area it is assumed that the variation 

of T (S ) is a linear function of r and z. The discretization will be carried out 
n n 

for the cosine modes only (T n (r,z,t)) the sine modes are very similar. 

For the purposes of this feasibility study the initial temperature was 
assumed uniform throughout the mirror, the sink for front surface radiation was 
uniform and the patches were heated one at a time. Within these restrictions the 
temperature distribution was symmetric with respect to a diametral plane through 
the center of the patch. Because of this the temperature and displacement fields 
were calculated with the patch at the correct radial position but centered about 
the x-axis. Then the temperature and displacement fields were rotated about the 
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r 


Figure 1 1-2 


z-axis until the patch was in the correct angular position. Then the displacement 
field was given a rigid body displacement to make the axial displacements zero 
at the supports. Because of this procedure it was only necessary to retain the 
cosine terms in the temperature field and the terms that are indicated later in 
the three displacement fields. 




(II-4) 


or 
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(11-5) 


T n = A ' “n 


an = AT V 


T = $ AT 


When I I -5 is substituted into I I -3a one has 
0 = 6T 


J $ 2 ^2 t ]^A AT n 


+ <ST t l\ t I pc$$ t rdA AT 
n I n 


{ 

I 


+ ST^A 1 J $p n rdS 


where ^ =(010) 

$ 2 * = (0 0 1 ) 

Identifying the integrals in I I -6 as C, D and q we have 


C n T n + D n T n + V" °’ 


(II -6 ) 


(1 1-7) 


The triangular elements are grouped by fours to form quadrilaterals as 
shown in Figure I I -3. 
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Figure 1 1 -3 

Here the numbers in triangles represent triangle corners, numbersin squares 
represent quadrilateral indicies, and the other numbers represent elements . 
The assembly of the finite elements are based upon two ideas: 

■\3i 

1. Global indexing of nodal point values related to local element nodal 
point values (Renaming) and 

2 . Summing the element nodal heat flux to equal the applied heat flux 
(Equilibrium) . 

The renaming is done by observation as follows 


T = T 1 = T 4 
1 1 2 


T = T 1 = T 2 
'2 '2 


The equilibrium requires 


q 2 = q 2 + q i 
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The reassembly is done via the renaming and equilibrium schemes operating on the sets 
of element equations 1 1 -7 . 

'After the quadrilateral elements have been forked from the triangular elements, 

the index 5 ; (middle index) is eliminated from the quadrilateral' representation. This 
is a straightforward elimination procedure and is not outlined here. 

The four-index quadrilateral elements are arranged to fill the r-z section of 
the mirror as indicated in Figure 1 1 -4 » 



Figure 1 1 -4 


Thus a node is located by I, J. In general the situation is as shown in Figure 1 1 -5. 
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Here 


q r j = q? + qj + q? + q!? 

= C 31 T I-1 , J-l + C 32 T I , 

+ C 4 1 T I , J-l + C 42 T I+1 , 
+ C ii T I , J + C 12 T I+1 , J 

. pb j + T 

L 2l'l-1 , J L 22 'i, J 
+ like DT terms . 

In this way the (1-1) x (J-l) quadrilateral 


j-l + C 33 T I, j + C 3 4 T I - 1 , J 

J-l + C 43 T I+1, J + C 44 T I s J 
+ C 13 T I+1 , J+l + C l*+T I ^ J+l 
+ C 23 T I , J+l + C 2 4 T I — 1 , J+l 

element equations of dimension 4x4 


are corresponding to each mode are compacted into one equation which is IM x JM in 
dimension. The coefficient matrices are banded matrices with a small number of 
elements on the band. 

The resulting equation 


CT + DT + q = 0 
n n n n M n 


(II -8) 


is integrated to determine the model temperature T r 
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Thermoelastic Formulation 

The temperature field T(r # e,z,t) determined from the previous thermal model 
as modal temperatures T n (r,z,t) are the -inputs, to the thermoelastic response. 

The thermoelastic response is the displacement field u- (x l5 x 2 , X 3 ,t), The 
inertial effects are ignored with the time dependence coming only from the time 
dependence of, the temperature field. The governing equations will be expressed 
in Cartesian tensor form although they are implemented in cylindrical coordinates. 
The displacements u^ must satisfy the Cauehy-Navier equations 


xu 


k,ki 


u(u. •• + -u . ■ . . ) - sT , . = 0 
M i,JJ *i 


(1 1-9) 


within the body and satisfy the boundary conditions 


!>k,k 6 ij + 


+ u j,i )]n j ’ 8Tn i 


( 11 - 10 ) 


on the surface of the body. Where X and y are Lame' elastic constants and @ is 
related to the linear coefficient of thermal expansion a by 


3 = a(3x + 2u) 


Equations ( 1 1 -9 and 10) are equivalent to the variational principle that potential 
energy must be stationary for equilibirum. This is given as 


■ <SPE 


J eTn ^ dS _ j eT^fiu. 

C 


dV 


't u k,k 6 ij + “‘ u i,j + 


dV = 0 


( 11 * 1 . 1 ) 
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The variational principle was expressed in' cylindrical coordinates and components 
of the displacement field u^, Ug, u z< This is not presented due to its length. 

The displacement fields were assumed as 

co 

u r (r,e,z) = 2 U n (r,z) Q0S n0 
n=0 

oo 

u 0 (r,0,z) = ^ V r ’ z) sin nQ 
n=l 

00 

u z (r,0,z) = ^ W n (r,z) cos -no 
n=0 

The sine terms for u f and u z and the cosine terms for u Q were not retained as 
previously explained. 

These displacement forms were inserted into the variational principle of 
Equation ( 11 - 11 ) to obtain individual variational principles for each mode (n). 

In addition the modal temperature and displacement fields are assumed to be linear 
functions of r and z inside each finite element triangle. Then, 


U 

V 

W 


= cx„ 


where = [U i Vi Wx U2 V 2 W 2 U3 V3 W3] 


and 


C = [Cil C 2 I C 3 I] 


c i = a, + b,r + d,z 
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Under these assumptions the variational principles become 


■ S PE n -0 = K( W .P n T n ) 


n = 0,1 ,, 


for an interior triangle, and 


' {PE n = 0 = Wn + P n T n ' R nV n = 0 * 1 - 


(11-13) 


(11-14) 


for a triangle having an edge on the surface. Where the coefficient matrices K^, 

P n , and R n can be identified from the integrals in the variational statement. 

The thermoelastic finite element equations for a single triangle are similar 
to the thermal equations. They are of dimension nine due to the three displacements 
of the three corners , however, there is no rate term. The assembly of triangles 
to form quadrilaterals and the systematic assembly of the quadrilateral elements 
proceeds exactly as in the case of the thermal program. For that reason it is not 
repeated here. The assembled equations are of the form 

KX + P T = 0 (H-15) 

n n n n 

Let ; M = IM x JM then K n is a square matrix of order 3M, X n is a row vector of order 

3M, P is 3M x M and T is a row vector of order M. 
n n 
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Response Computation 

The thermal and deflection programs whose -modal -finite element discretization 
have been discussed in the two previous sections yield the following vector matrix 
equations 

CT + DT = q for each mode 

and 

KX + PT = 0 for each mode 

where 

q is the vector of nodal heat inputs 
T is the vector of nodal temperature 
X is the vector of nodal displacements „ 

To compute the displacement of any node at time t the thermal equation is 
integrated to yield the temperature and the displacement is computed from the 
displacement equation where the nodal heat input rates q must be inputed. 

The integration algorithm that was used is as -follows 



where h is the integration delta time. The thermal equation via the integration 
algorithm yields 


(C + 


2Dw T i+l H 
h M 2 


T. 

i 


l DT i 


2<di 


( i+l 


i+1 


2 ( 


i+1 


+ T. 


- 1 ) - 


T. 

l 
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In using this procedure to iterate the nodal temperatures , no problems in 
matrix inversion were encountered* The anticipation q-j+i was removed by assuming 

■ vta + Vi)’ 

a reasonable assumption for slowly changing heat -inputs and high iteration rates. 

The displacement vector of the body is obtained from the nodal temperature 
vector by matrix inversion and multiplication in the thermo-elastic equation. No 
problems with matrix inversion were encountered. 

When the deflection associated with each successive mode is computed it is 
added to the preceeding modest The vector displacement of each node is established. 
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The Influence Matrix Computation 

The influence matrix computation is associated with the complete steady T state 
response of the nodes of .the body excited by a designated pattern of constant. heat 
inputs . A subset of the nodes on the front surface of 'the mirror is assigned the 
role of output nodes. A pattern of patches is specified on the rear of the mirror. 
A unit heat rate is supplied to the patches one at a time and the vector of 
longitudinal steady state deflection is recorded. These deflection vectors are 
the columns of the influence matrix A, where 


w = Aq 

w is the longitudinal output deflection vector and q is the heat-input vector 
whose components represent the heat input at each patch. 

This influence matrix may be.used in control simulation for any output node 
set and patch pattern set which are subsets of the set for which the influence 
matrix is computed. 
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Removal of Free Thermal Expansion 

In heating the mirror for control purposes a significant amount' of heat goes 
into the free thermal expansion of the' mirror. This deflection does not contribute 
to local error smoothing and may be compensated by refocussing of the mirror. 
Moreover, if retained as part of the control, it significantly reduces the 
sensitivity of the control function. For these reasons the free thermal expansion 
terms are measured and their effect is removed from the influence matrix. 

Consider the mirror geometry shown in Figure 1 1 -6 . 
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The difference in the two reference spheres is 


y R > 


The best fit sphere is found by minimizing the functional 


J = (w. ... - y . 

ij J 


where w. . is the displacement of the ij— surface nodal point due to heat input 
J 

on the rear. 

The minimization yields 


‘ = "r R 


cos 0,. 

•J w ijfcos.-0“ " ^ 

COS 0. ~~ 

M - 1]2 
1 0 


which represents the adjustment of the reference sphere, 
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Error Computation and Control 

The deflection of the front of the mirror is composed of a disturbance w Q and 
control induced deflection 

W c = Aq 


so 


w = w Q + Aq 

The front surface error figure is 

J - w^Qw , 

where Q is a symmetric weighting matrix. Since one may wish to minimize the use 
of the control a control effort cost term 

C = q t Rq 

is also introduced, where in most cases R is a diagonal matrix, whose elements 
reflect the priority associated with control cost versus surface error reduction. 
The minimization of the index 

w^Qw + q^Rq 


under the constraint 


w - WpAq 


yields the optimal control vector 

q = =[A t QA'+ R ]" 1 A t Qw D 
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and associated best index 



III. DISCRETION OF FORTRAN PROGRAMS 
A. RESPONSE 


Purpose 

This program computes the thermoelastic deflection of a set of nodes specified 
in an axisymmetric body. Heat is inputed on one surface of the body and radiates 
at the other. 

Options 

This program provides the following options. 

(1) The temperature and deflection of the node set may be computed at any time. 

(2) The steady state temperature and deflection may be computed. 

(3) The influence matrix of surface deflection for a designated heat input; 
pattern may be computed. 

Input Parameter Definition 
Parameter 
NUMPAT 
IM 

JM 

CK 
CP 

CF 
TO 


Definition 

The total number data cases to be run. 

Number of node locations in the radial direction. Cannot 
exceed 15. 

Number of node locations in the logitudinal direction, 
Cannot exceed 5. 

Thermal conductivity of the material. 

Specific heat of the material (0. value indicates steady- 
state solution) . 

Radiation coefficient of the front surface. 

Temperature of radiating reference medium. 
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Parameter 


Definition 


DELT 

NTS 

IPRINT 

TI 

ALF 

E 

V 

NM 

IS 

KM 

52 

53 

INFLU 


IX 


IP 

KP 

PA 


Integration time. 

Number of integration steps. 

Number of integration steps between printings. 

Initial mirror temperature. 

Thermal expansion coefficient, 

Young’s modulus. 

Poission's ratio.. 

Highest harmonic (angular) component. 

Location of radial node of support ring. 

Number of angular divisions. Must not exceed 12. 
Angular position of second support (first ope is at 0) . 
Angular position of third support. 

Switch to choose simulation of response or influence 
matrix computation. 

If INFLU = 0 Program calculates the response to heating 
one patch. 

If INFLU > 0 Program calculates the thermoelastic in- 
fluence coefficient matrix. This matrix 
is written in a data file through 1-0 
unit 4. 

Index of the radial ring on the back at which the modal 
boundary conditions are applied. 

Radial node of single patch input. 

Angular position of single patch input. 

Patch angle for single patch. 
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Parameter 

PH 

PAN (I) 

DO 

DI 

H 

FNO 


Definition 

Heat input rate for single patch. (Negative is input 
positive is output) , 

Patch angle for i-th radial node in full pattern 
generation . 

Outside diameter of mirror. 

Inside diameter of mirror. 

Thickness of mirror. 

F-number of the mirror = Focal Length/ Diameter of Mirror. 
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Input Data Card Listing 


Card No. 

Parameter 

Data Field 

Format 

1 

NUMPAT 

1-5 

15 

2 

IM 

1-5 

15 


JM 

6-10 

15 


CK 

11-20 

F10.5' 


CP 

21-30 

F10.5 


CF 

31-40 

F10.5 


TO 

41-50 

F10.5 


DELT 

51-60 

F10.5 


NTS 

61-65 

15 


IPRINT 

66-70 

15 


TI 

71-80 

F10.5 

3 

ALF 

1-10 

F10.5 


E 

11-20 

F10.5 


V 

21-30 

F10.5 

4 

NM 

1-5 

15 


IS 

6-10 

15 


KM 

11-15 

15 


S2 

16-20 

15 


S3 

21-25 

15 


INFLU 

26-30 

15 


IX 

31-35 

15 
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Card No . 

Parameter 

Data Field 

Format 

5* 

IP 

1-5 

15 


KP 

6-10 

15 


PA 

11-20 

F10 


PH 

21-30 

F10 

5** 

PAN(l) . . .PAN (7) 

1-70 

7F10.5 


5+N 

PAN(IM-7) . . .PAN(IM-l) 

1-70 

7F10.5 

6+N 

DO 

1-10 

F10.5 


DI 

11-20 

F10.5 


H 

21-30 

F10.5 


FNO 

31-40 

F10.5 


Cards 2 through 6+N constitute one data set 
and there should be NUMPAT data sets. 


*,■ Used only if INFLU equals zero. 

** Used only if INFLU is greater than zero. 
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Output of Program 

A. Option INFLU = 0 

1. Repeated input data. 

2. Table of temperatures and deflection of the node (I=IM, J=JM) for each mode 
for a check .of convergence. 

3. Table of displacements and temperatures of the nodes on the front surface 
before the support constraints are added. 

4. Table of displacements of the nodes on the front surface after supports are 

added. 

B. Option INFLU = 1 

1. Repeated input data. 

2. Patch angles of each patch ring. 

3. The displacements of the nodal points on the front surface of the mirror 
corresponding to each patch location are read into a file. 
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Heat Input and Deflection Measurement-Points . 

The node distribution for both the front and the rear of the mirror is shown 
in Figure III-l. The surface has IM concentric nodal rings and KM nodal rays. 

The nodes occur at intersections of rings and rays. The nodes are selected by the 
program after mirror dimensions, KM, and IM are given. 

The nodes for measurement of surface deflection are chosen in the control 
program. Any subset of the full set of nodes may be chosen. 

The heat is applied to the rear by patches. The patches extend the entire 
distance between rings and make an angle PAN centered about a ray. Patch location 
and patch angle may be inputed in the control program. 

Several typical patch location patterns are shown in Figure IIIr2 through 
II 1-5. 
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Sample Points on Front Surface 


Figure III-l 
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Heat Patch Locations 
-Typical Patch- 


Figure 1 1 1-2 
-32r 










12 Patch Distribution 


Figure II 1-5 
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B. CONTROLS 


Purpose 

This program computes the patch heats necessary to minimize the weighted surface, 
error of the mirror taken at designated sample points, using designated heater 
locations. In addition it computes the surface error both before and after applica- 
tion of the thermal input. The disturbance error is computed internally. 


Input Parameter Definl tion 

Parameter 

IM 

KM 

DI 

DO 

FNO 

NHP 

NSP 

IHP 

ISP 

A(I,J) 


Definition 

Number of nodes in raidal direction 

Number of angular divisions 

Inner Diameter 

Outer Diameter 

Focal Length/Diameter 

Number of heater locations 

Number of sample points 

Individual heater location 

Individual sample point 

Coefficients of the influence matrix 

computed in RESPONSE (read from file) 
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Card No. 

Parameter 

Data Field 

• Format 

1 

• IM 

1-5 

15 

1 

KM 

6-10 

15 

1 

DI 

11-20 

F10. 

1 

DO 

21-30 

F10. 

1 

FNO 

31-40 

F10. 

2 

NHP 

1-5 

15 

2 

NSP 

6-10 

15 

3 

IHP(I) 

1-80 

1615 

4 

ISP(I) 

1-80 

1615 


Output of Program 


1. Repeated heat patch points, 

2. Repeated sample points, 

3. Coefficients of the reduced influence matrix corresponding to the 
heat patches and sample points selected. 

4. The surface error of the sample points before control. 

5. The performance index before control . 

6. The performance index after control . 
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IVv FEASIBILITY RESULTS 


The results presented here were obtained for a small fused-silica mirror 
whose properties and dimensions are given in Table IV-1. The deflection patterns 
are shown in Figures IV-1 through IV-10 for various one patch heater locations. 
The mirror is supported at 120 5 ® intervals as shown in the figures. Figures IV-1 
through IV-6 show the steady state deflection patterns for heat inputs of 3 watts 
and various patch locations. Figures IV-8 and IV— 9 shows a, transient deflection 
pattern for the same heat rate. Figures IV-10 and IV-1 1 show the transient de- 
flection of selected points on the front surface as a function of time. These 
points are shown in Figure IV-7. 
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SMALL FUSED SILICA MIRROR 


“INNER- RADIUS 0. 

'OUTER RADIUS 20. IN 


•THICKNESS 


.667 IN 


Structural Properties 


•SPECIFIC HEAT 
'CONDUCTIVITY 
•EMISSIVITY 
•POISSON'S RATIO 
’YOUNG'S MODULUS 

“COEFFICIENT OF THERMAL EXPANSION 


.015 BTU/IN 3 DEG F 
o 05 BTU/HR*IN* DEG F 
.04 
.17 

106 10 s LB/ IN 2 
.311*10~ 6 /DEG F. 


Response- Properties 

- ’NOMINAL HEATING POWER 3 WATTS 

•THERMAL TIME- CONSTANT 19 HOURS 


Table IV-1 
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Transient Response 
Heat Input 10" 3 BTU/HR'IN 2 

Figure IV-7 
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ZERO DEFLECTION CONTOUR 
ONE HOUR 

Figure IV-8(a) 
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10 


ZERO DEFLECTION CONTOUR 
TEN HOURS 
Figure IV-8(b) 
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ZERO DEFLECTION CONTOUR 
FIFTY HOURS 


Figure IV-8(d) 



1 y MICRO- INCH CONTOUR 
FIVE HOURS 

Figure I V-9 ( a) 
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1 y INCH CONTOUR 


TEN HOURS 
Figure I V-9 ( b) 
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- APPENDIX 


PROGRAM LISTING 
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o o o ooooooooooo o o ooooooooo 


main 

PROGRAM TO DETERMINE TEMPERATURE DISTRIBUTION AND 

thermal distortion of a mirror 
parameters 

CP=SPECIFIC HEAT 
CK=THERMAL CONDUCTIVITY 
CF=RADlATION FACTOR FOR FRONT FACE 
TO=REFERENCE TEMPERATURE 
I»J=ROW AND COL INDICES OF GRID PT 
R(I» J) #Z(I# J)=COORDlNATES OF GRID POINT 
IM»JM=MAX VALUES OF I AND J 
NC=NUMBER OF THE HARMONIC COMPONENT 

T C ( I ) =temperature component 

DELT= INTEGRATION TIME STEP 

NTS=NUM8ER OF TIME INTERVALS 

ALF=LINEAR COEFFICIENT OF THERMAL EXPANSION 

E=YOUNGS MODULUS OF elasticity 

V= POISSONS RATIO 

D0=0UTSIDE DIAMETER OF MIRROR 

DI-INSIDE DIAMETER OF MIRROR 

HrTHlCKNESS OF MIRROR 

FNO=F-NUMBER OF THE MIRROR (FOCAL LENGTH/DO) 

DIMENSION AS ( 15# 15) »DK(15#15) »DKK(15» 15) »DKKK(15»15) » IPIV0T(15) 
lU( 15# 15) » INDEX (15# 2) #ZZ(15»1) »SPP(15) »SPPP(15) #Q2(15»15) # 
2SAP(15) »SS(15) » SSS ( 15) »SP(15#15) #P1 ( 15# 15# 15) 

DIMENSION D ( 75) »R(15»5) #Z(15#5> »C (75#75#2) rQR(15) #QB(15#14) , 
lTC ( 75) #X(75) r A (75 #75) #Q(75) » 01(75) »P(14) »F(15»15) 

DIMENSION WF(15#12) »WB(l2) #TT(15»12) 

DIMENSION HEAT ( 100 ) 

DIMENSION PAN(14) 

INTEGER SI# S2» S3 
COMMON R»Z 

COMMON/STIF/KD # KS # PD » PS » PL # A 
COMMON /CND/C#D»QR»QB 

DIMENSION KD(15»15»15) »KS(15» 15# 15) # PD (15# 5 #15) # PS (15# 5 #14 ) , 
3PL(15»5»14) »KL(15»15»15) 

REAL KD#KS»KL 

LOOP TO INPUT ALL THE PATCHES 

READ (5»3100)NUMPAT 
3100 FORMAT (15) 

DO 72 NU=1#NUMPAT 
WRITE(6#6500) 
t>500 FORMAT (1H1) 

READ(5#iOO)lM»JM»CK#CP»CF#TO»DELT»NTS#IPRINT#TI 

100 FORMAT (2I5#5F10»5#2I5#F10«5) 

READ ( 5# 101 ) ALF »E # V 

101 FORMAT (3F10«5) 

WRITE (6# 102) IM» JM»CK»CP»CF»TO#DELT»NTS 
WRITE(6» 107) ALF»E» V 
READ (5# 2600 >NM# IS#KM»S2#S3» INFLU# IX 
2600 F0RMAT(7I5) 

nmp=nm+i 

WRITE (6 #2800 ) TI 

2800 format (27H initial temperature# ti = »fio.s/////) 
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WRITE (6 > 2601 ) NMP>KM> IS>S2>S3>IX 

2601 FORMAT (33H NUMBER OF HARMONIC COMPONENTS = >15/ 

137H NUMBER OF ANGULAR PATCH POSITIONS = >15/ 

231H RADIAL NODE OF SUPPORT RING = >15/ 

348H ANGULAR POSITIONS OF SUPPORTS ARE>K =. 1* K = >I2>2X> 

410H AND> K = >12/ 

56H IX = >15/) 

Nl=3 

MM=15 

LL=15 

MAT=IM 

fkm=km 

NTP=NTS+1 

IM1=IM-1 

M=IM*JM 

FJM1=JM-1 

FIM1=IM1 

IF(INFLU.GT.O)GO TO 2500 
READ <5» 2602) IP >KP>PA>PH 

2602 FORMAT(2I5>2F10.5) 

WRITE (6> 2603) IP>KP>PA> PH 

2603 FORMAT (52H THERMOELASTIC RESPONSE OF THE MIRROR IS CALCULATED/ 
124H RADIAL NODE OF PATCH = >15/ 

229H ANGULAR POSITION OF PATCH = >15/ 

315H PATCH ANGLE = >F10.5/ 

414H PATCH HEAT = >F10.5) 

GO TO 2501 

2500 READ<5>2604) (PAN( I ) > 1=1 > IM1 ) 

2604 FORMAT (7F10 • 5) 

WRITE (6>2605) (I>PAN(I) >1=1 >IM1) 

2605 FORMAT (72H THERMOELASTIC INFLUENCE COEFFICIENTS ARE CALCULATED AND 
1 WRITTEN ON FILE/ 

24X>1HI>5X>11HPATCH ANGLE/ 

3(I5>5X>F10.5) ) 

2501 CONTINUE 
PI=3. 1415927 

GRID GENERATION FOR SPHERICAL MIRROR 

READ (5> 112) DO>DI>H>FNO 
WRITE (o>113) DO>DI>H>FNO 
IF(FN0.GT.99) GO TO 201 
RF=2.*D0*FN0 
RB=RF+H 
DR=H/FJM1 
SI=.5*DI/RF 
CI=SQRTll.-SI*SI) 

S0=.5*D0/RF 
CO=SQRT ( 1 . -SO*SO ) 

THETI=ATAN(SI/CI) 

THETO=ATAN (SO/CO) 

DTHET=( THETO-THETI ) /FIMl 
EM=THETI 
DO 39 1=1 > IM 
RR=RB 

DO 38 J=1 > JM 
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R(I»J)=KR*SIN(EM) 

. Z(I»J)=rs-RR*COS(EM) 

36 RR=RR-Dk 
39 EM=EM+DTHET 
60 TO 202 

201 DH=.5*(00-DI)/FIM1 
DZ=H/FJM1 
DO 200 1=1 » IM 
FI=I-1 

DO 200 J=1»JM 
FJ=J-1 

R ( I * J) =0 . 5*DI+FI*DH 
200 Z(I»J)=FJ*DZ 
202 CONTINUE 
C 

c call subroutine to calculate c»d»kd>ks»pd»ps»pl 
c 

CALL CAND (iMr JM»M»CK»CPrCF) 

c 

components of patch heat 

c 

IF l INFLU • LE • 0 ) GO TO 250 2 
IP=1 

2505 PA=PAN(1P) 

IP1=IP+1 

AP=PA*3 . 14159265/360 . * ( R ( IP1 » 1 ) **2-R ( IP * 1 ) **2 ) 

PH=-1./AP 
2502 CONTINUE 

DO 2000 N=1»NM 
FN=N 

HEATO=PH*PA/360. 

2000 HEAT(N)=C2.*PH*SIN(FN*Pa*PI/360.) )/(FN*PI) 

C 

C INITIALIZE WF AND WB 
C 

DO 2001 K=1»KM 
WB(K)=0.0 
DO 2001 1=1 » IM 
TT ( I » K ) =0 • 0 

2001 WF ( I »K) =0 • 0 
C 

C START MODE LOOP 
C 

IF(INFLU.GT.O) GO TO 6000 
WRIT£(6» 500) 

500 FORMAT (/////31H CONVERGENCE OF MODAL RESPONSES// 

11X»4HM00E» 11X» 10HHEAT INPUT r 4X» 1HI r4X» 1HJ, lOXr 10HDEFLECTION, 
29X » HHTEMPERATURE ) 

6000 DO 2002 NT=1»NMP 
Nl=3 
NC=NT-1 

IF (NC • EQ « 0 ) Nl=2 
JM3=N1*JM 
LCV=JM3 
I JM=lM*JM3 

CALL STIFF (E » V » ALF r IM» J m»NC » CK »CF) 
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non 


REARRANGE KD AND KS AND FORM KL FOR POTTER 

DO 50 I=1»IM 
DO 50 Ll=l#JM3 
DO 50 L2=1»JM3 

50 KL(L1»L2»I)=KD(L1»L2»I> 

DO 51 I=1»IM 

DO 51 Li=l»JM3 
DO 51 L2=1»JM3 

51 KD ( I » Ll » L2) =KL (LI » L2» I ) 

DO 52 1=1 »IM 

DO 52 Ll=l»JM3 
DO 52 L2=1»JM3 

52 KL (LI »L2» I ) =KS (LI »L2» I ) 

DO 53 1=1 » IM 

DO 53 Ll=l t JM3 
DO 53 L2=l t JM3 
53 KS( I 'Ll »L2)=KL(L1 »L2» I ) 

DO 54 I=2» IM 
DO 54 Ll=l » JM3 
DO 54 L2=1»JM3 

54 KL(I»Ll»L2)=KS(I-l»L2*Ll> 

c initial conditions on temperature 

c 

IF(NC.NE.O) GO TO 131 
DO 37 1 = 1 » M 
37 TC ( I ) =TI 
GO TO 132 

131 DO 133 1=1 » M 
133 TC(D=0.0 

132 CONTINUE 
C 

C INITIALIZE DISPLACEMENTS 

C 

DO 24 L1=1»IM 
DO 24 L2=1»JM3 
24 U(L1'L2)=0.0 
C 

C RADIATION FROM FRONT 

C 

FC=NC 

DO 20 Ll=l t M 
X (Ll ) =0 • 0 
20 Q1(L1)=0.0 

IF(NC.NE.O) GO TO 22 
DO 21 Ll=l » IM 
LR=L1*JM 

21 Q1(LR)=Q1(LR>+ T0*QR(L1) 

22 CONTINUE 
C 

C SET UP AND INVERT EFFECTIVE CONDUCTIVE MATRIX IN A 
C 

DO 2 Ll=l » M 
DO 2 L2=l »M 

A(L1»L2)=C(L1»L2»1)+FC*FC*C(L1»L2»2) 
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IF(L1.EQ.L2)A<LDL2)=A(LDL2)+2.0*D(L1)/DELT 
2 CONTINUE 

CALL INVERT (A»Mf 75) 

INITIALIZE SOLUTION AND HEAT INPUT 


DO 99 1 = 1 » IM1 
99 P(I)=0.0 


STEP BY STEP SOLUTION 
IPR=0 

IF(CP.6T. 0.0001) GO TO 73 

NTP=2 

IPRINT=1 

73 DO 6 I 1=1 »NTP 
DO 45 Kl=l»M 
45 X (K1 ) =0 • 0 

FI=I 1-1 
TIME=Fl*DELT 
IFdI.EQ.DGO TO 40 
IPR=IPR+1 

HEAT INPUT ON BACK 

IF(NC .NE. 0) P(IP)=HEAT(NC) 

IF (NC .EQ. 0) P(IP)= HEATO 
DO 26 L=2» IM1 
I J= (L-l ) *JM+1 

28 Ql(IJ)=QB(L»L-l)*P(L-l)+QB(LfL)*P(L) 
Q1 ( 1 ) =QB (1»1)*P(1) 

I2=M+1~JM 

Q1(I2)=GB(IM»IM1)*P(IM1) 

SET UP THE COMBINE HEAT INPUT VECTOR 
DO 3 Ll=l »M 

3 Q(L1)=2.0*D(L1)*TC(L1)/DELT +Q1(L1) 

CALCULATE THE NEW TEMPERATURE 

o002 DO 10 Li=l t M 
DO 10 L2=l t M 

10 X(L1)=X(L1)+A(L1»L2)*Q(L2) 

IF(CP.LT. 0.0001) GO TO 70 
DO 4 Ll=l»M 

4 TC(L1)=-TC(L1)+2.*X(L1) 

IF ( IPR.NE. IPRINT) GO TO 6 
IPR=0 

5001 CONTINUE 

IF (CP. GT .0.0001) GO TO 5000 
70 DO 62 Ll=l t M 
62 TC(L1)=X(L1) 


C 
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SET UP THE THERMAL FORCE VECTOR 


C 
C 

aOOO IF(NC.NE.O) GO TO 71 
DO 55 Li = l,M 
55 TC(L1)=TC(L1)-T1 
71 DO 41 Ll=l , IM 
DO 41 L2=1,«JM3 

41 F(L1»L2)=0.0 
DO 43 1 1=1 , IM 
DO 43 L2=1»JM 
DO 43 J=1 , JM 
DO 42 K=1,N1 
Ll=Nl*(J-l)+K 
L3=JM*(il-l)+L2 
L4=L3+JM 
L5-L3-JM 

F(I1,L1)=-PD(L1,L2,I1)*TC(L3>+F(I1,L1> 
IF(I1.N£.IM>F(I1,L1>=F(U»L1)-PS(L1»L2,I1)*TC<L4) 

IFdl.EQ.DGO TO 42 

F(I1,L1)=F(I1,L1)-PL<L1,L2»I1-1)*TC(L5) 

42 CONTINUE 
43 CONTINUE 

IMPOSE BOUNDARY CONDITIONS 
N=NC 

IF (N»NE. 0 ) GO TO 49 
DO 44 Ll=l,JM3 
IF(IX.NE.1)KS<IX-1»L1»2)=0. 

IF<IX.NE.IM)KL(IX+1,L1,2)=0. 

KD ( IX , Ll » 2 ) =0 • 

KS ( IX » 2 » Ll ) =0 • 0 
IF(IX.NE.1)KL(IX»2»L1)=0.0 
44 KD(IX»2rLl)=0.0 
KD<IX,2,2)=1.0 
F ( IX »2) -0 • 0 
GO TO 48 

49 IF (N»NE • 1 ) GO TO 48 
DO 47 L2=l»3 
DO 46 L1=1»JM3 
KS(IX,L2,L1)=0.0 
IF(IX.NE.1)KL(IX,L2»L1)=0.0 

46 KD(IX»L2»L1)=0.0 
KD(IX,L2,L2)=1.0 

47 F(IX»L2)=0.0 

48 CONTINUE 
DO 9000 L1=1,JM3 
DO 9000 1.2=1 »JM3 

9000 AS(L1,L2)=KL(2,L1»L2) 

WRITE (6, 9002) ( <F(L1,L2) ,L2=1,15) ,L1=1,15> 

9002 FORMAT (15E8* 3) 

CALL MATINV(AS,JM3,SS»0,DETM,IPIVOT,INDEX,MM,ISCALE> 

WRITE (6,6001 )N»DETM 

oooi format (5H n = »I2,5x,7hdetm = ,E20.8) 
calculate displacements 

call POTTER (KL, KD,KS,F,U, MAT, LCV, AS, DK,DKK,DKKK,IPIVOT, INDEX, ZZ, 
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4SPP » SPPP f Q2 1 PI » SAP » SS t SP t SSS t MM » LL ) 

40 CONTlNUc. 

WRITE (6 >9003) < <U(L1>L2) >L2=1»15) »L1=1»15) 

9005 FORMAT ( 7H U( I > J) / ( 15E8.3) ) 

6 CONTINUE 
C 

COMPUTE WB AND WF 
C 

IF (NC .NE. 0 ) GO TO 4000 
DO 4001 L1=1>M 
4001 TCCL1)=TC(L1)+TI 
4000 CONTINUE 

DO 2003 K=1>KM 

FN=NC 

FK— K 

WB(K)=WB(K)+U(ISrND*C0S(FN*(FK-l.)^PI*2./FKM) 

DO 2003 1=1 > IM 
IT=I*JM 

TT ( 1 > K ) =TT ( I r K ) +TC ( IT ) *COS ( FN* ( FK-1 . ) *PI*2 . /FKM ) 
2003 WF(I>K)=WF(I>K)+U(I t JM3)*C0S(FN*(FK-1. )*PI*2./FKM) 
C 

C OUTPUT RESULTS OF A MODE 
C 

IF (INFLU. GT.OJGO TO 2002 
1 = 1 
J=JM 
JJ=J*Nl 

JJJ=JM*(I-1)+J 

WRITE (6# 501) NC>P(IP) *I> J>U(I> JJ) >TC(JJJ> 

501 FORMAT ( I5>E20«8>2I5>2E20»8) 

2702 CONTINUE 
2002 CONTINUE 

DTHD=360./FKM 
IF(INFLU.GT.O)GO TO 3000 
WRITE (6> 1501) 

DO 1006 K=1 > KM 
WRITE (6> 1505) 

FK=K-1 

thd=fk*dthd 

DO 1006 1=1 > IM 

1006 WRITE (6# 1502) I>K>R(I»JM) »THD>WF(I>K) »TT(I»K) 

WRITE (6> 1503) 

DO 1007 K=1 >KM 
FK=K-1 

thd=fk*dthd 

1007 WRITE (6 > 1504) K >RS> THD> WB IK) 

C 

C ROTATION OF OUTPUT FOR SUPPORTS 
C 

5000 CONTINUE 

DIMENSION WBR(12)»WFR(15»12) 

Sl=l 

1F(INFLU.LE.O)GO TO 2503 
KP=1 

2503 CONTINUE 
FKM=KM 
FS2=S2 
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FS3=S3 

DTHET=2.*PI/FKM 

DTHD=36U./FKM 

KR=KP-1 

DO 1001 K=1»KM 
KR=KR+1 

IF (KR.EQ. KM+1 ) KR=1 
DO 1002 1=1 » IM 

1002 WFR(I»KR)=WF(I»K) 

1001 WBR(KR)=WB(K) 

1000 CONTINUE 

W1=WBR(S1) 

W2=WBR(S2) 

W3=WBR (S3) 

RS=R(IS»1) 

THET2= (FS2-1 . ) *DTHET 
THET3= (FS3-1 . ) *DTHET 
X1=RS 

X2=RS*C0S(THET2) 

X3=RS*COS(THET3) 

Y2=RS*SIN(THET2) 

Y 3=RS*SIN ( THET3 ) 

DET=Y3*(X2-X1)+Y2*(X1-X3) 

IFtDET.bT.O.OOODGO TO 1003 
WRITE (6# 1510 ) 

1510 FORMAT (23H SUPPORTS LOCATED WRONG) 

GO TO 72 

1003 WO=( (X2*Y3-X3*Y2)*W1-X1#Y3*W2+X1*Y2*W3)/DET 
THX= ( ( X3-X2 ) *W1 + ( X1-X3 ) *W2+ ( X2-X1 ) *W3 ) /DET 
THY=( (Y3-Y2)*W1-Y3*W2+Y2*W3)/DET 

IF (INFLU. GT.O) GO TO 3001 
WRITE (6# 1500 ) WO»THX»THY 

3001 CONTINUE 

1500 format (//35H rigid body displacement parameters/17h z-translation 

1= »E15.B»5X>13HX-R0TATI0N = r E15.8* 5X r 13HY-R0TATI0N = *E15.8//) 

DO 1004 K=1»KM 
FK— K 

DO 1005 1=1 * IM 
RI=R(I» JM) 

TK=(FK-1.0)*DTHET 

XIK=RI*C0S(TK) 

YIK=RI*SIN(TK) 

1005 WFR(I»K)=WFR(I»K)-W0-YIK*THX+XIK*THY 
XIS=RS*COS(TK) 

YIS=RS*SIN(TK) 

1004 WBR(K)=WBR(K)-WO-YIS*THX+XIS*THY 
IF ( INFLU »GT . 0 ) GO TO 3002 

WRITE (6*1501) 

DO 1506 K=1 »KM 
WRITE(6* 1505) 

FK=K-1 

thd=fk*othd 

DO 1506 1=1 » IM 

1506 WRITE(6»1502)I»K»R(I»JM) »THD*WFR(I rK) *TT(I*K) 

3002 CONTINUE 

IF ( INFLU »LE • 0 ) GO TO 72 
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WRITE(6*5Q99) IP*KP 
o099 FORMAT (6 H IP = *I5*5X»5hKP = *15) 

WRITE (4 » 1600) ( ( WFR ( I » K ) , 1 = 1 * IM) * K=1 »KM) 
1600 FORMAT (bE13*8) 

WRITE (4 *7006) ( (WFR(I*K) *I=1»IM> »K=1*KM) 
7006 FORMAT (bE13*8) 

7005 CONTINUE 

IF (INFLU. GT.O) GO TO 3003 
WRITE (b» 1503) 

DO 1507 K-l » KM 


FK=K-1 

THD=FK*DTHD 

1507 WRITE(6,1504)K*RS*THD»WBR<K) 

1501 format 129 H z-displacements on the front //4X»ihi»4x»ihk»iox» 

1 1HR » 14X * 5HTHET A » 12X * 1HW » 16X » 11HTEMPERATURE) 

1502 FORMAT(2I5»2E15.5*2E20.8> 

1503 format (//5oh z-displacements at the support radius on the back// 

l4X * 1HK » 14X* 2HRS# 13X* 5HTHETA * 10X# 1HW) 

1504 FORMAT ( 15* 3E20 *8) 

102 F0RMAT(////39H TEMPERATURES AND THERMAL DISPLACEMENTS*/// 

119H NO* OF GRID ROWS = *I5*5X»20H NO* OF GRID COLS* = *15/ 

2/18H MIRROR PROPERTIES/ . 

320H CONDUCTIVITY r CK = ,E15*6/21H SPECIFIC HEAT* CP = »E15*6/ 

424H RADIATION FACTOR* CF = *E15.6/29H REFERENCE TEMPERATURE* TO = 
5»E15*6//24H SOLUTION TIMER C0NTR0LS/32H INTEGRATION TIME STEP* DELT 
bTA = *E15.6*5X*28HNUMBER OF TIME STEPS* NTS = *15//) 

107 FORMAT (48H LINEAR COEFFICIENT OF THERMAL EXPANSION* ALF = »E15*6/ 
135H YOUNGS MODULUS OF ELASTICITY* E = *E15*6/ 

321H POISSONS RATIO* V = *E15.6/////) 


112 FORMAT (4F10 • 5) _ # 

113 FORMAT (7///16H MIRROR GEOMETRY/ /20H OUTSIDE DIAMETER - »F10.5/ 

119H INSIDE DIAMETER = »F10*5/13H THICKNESS = *F10.5/12H F-NUMBER = 
2 »F10.5///) 

IF ( INFLU.LE* 0 ) GO TO 72 
KP=KP+1 

IF (KP.GT .KM) GO TO 2504 
GO TO 2503 
*504 IP=IP+1 

IF ( IP.LE. IM1 ) GO TO 2505 
72 CONTINUE 
STOP 
END 
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CAND 

SUBROUTINE CAND( IM» JM»M,CK»CP»CF) 

COMMON R»Z 

COMMON /CND/C»D»QR>QB 

DIMENSION C(75f 75#2) »D(75) »R(15»5) »Z(15»5) 

DIMENSION CT(3»3»2)»CQ(5»5»2)»6T(3»3»;2)*DQ(4)»A(3*3)»F1(5»5*2)» 
2F2(5»5»2) »G(5»5#2) »D1(5) »D2(5) »GQ(3) »QR(15) »QB(15»14) 

DIMENSION RZ(3> »ZR(3) rX(10) 

IM1=IM-1 
JM1=JM“1 
DO 20 Ll=l>M 
D(L1)=0.0 
DO 20 L2=1#M 
DO 20 L3=l»2 
20 C(L1»L2»L3)=0.0 
DO 47 L=1»IM 
47 QR(L)=0.0 

DO 60 M1=1»IM 
DO 60 M2=1»IM1 
60 QB (Ml » M2) = 0 • 0 

DO 25 L1=1»JM 
Dl(Ll)=Q.O 
D2(L1)=0.0 
DO 25 L2=l t JM 
DO 25 L3=l»2 
F1(L1»L2»L3)=0.0 
F2(LlrL2»L3)=0.0 
25 G(L1»L2»L3)=0.0 
DO 24 1=1 »IM 
DO 46 L=1 t JM 
D2 (L) =0 . 0 
DO 46 LL=1 » 2 
DO 46 K=1 t JM 
G(L»K»LL)=0.0 
46 F2(L'K>LL)=0.0 
DO 21 J=1 » JM1 
IF (I.EQ.IM)GO TO 21 
R1=R ( I » J) 

Z1=Z(I»J) 

R2=R<I+1» J) 

Z2=Z ( I + l » J) 

R3=R(I+1» J+l) 

Z3=Z(I+Jl»J+1) 

R4=R(I»J+1) 

Z4=Z(I»J+1) 

R21=R2-R1 

R32=R3-R2 

R41=R4-R1 

R34=R3-R4 

Z21=Z2-Z1 

Z32=Z3-Z2 

Z41=Z4-Z1 

Z34=Z3-Z4 

AREA=.5*(R41*Z41-R21*Z2l-R32*Z32+R34*Z34)-R32*Z21+R34*Z41 
AR=.5*(R41*Z4l)*(Rl+2.*R41/3.0)-.5*R21*Z2l*(Rl+2.*R21/3. )-R32*Z21* 
3 ( R2+ • 5*R32 ) - . 5*R32*Z32* ( R2+2 . *R32/3 • ) + • 5*R34*Z34* ( R4+2 * ♦R34/3 • ) +R3 
44*Z4l*(R4+.5*R34) 
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AZ=« 5*R41*Z41* (Z1+Z41/3. )-»5*R21#Z21* ( Z1+Z21/3. ) -R32*Z21* (Zl+«5*Z2 
bl ) 5*R32*Z32* ( Z2+Z32/3 . ) ♦ . 5*R34*Z34* ( Z4+Z34/3 . ) +R34*Z41* ( Zl + . 5*Z4 
bl) 

RC=AR/Ar(£A 
ZC=AZ/AREA 
DO 16 Ll=l t b 
DO 16 L2=l t b 
DO 16 L3=l»2 

16 CQ(L1#L2»L3)=0.0 
DO 17 Ll=l t 4 

17 DQ(L1)=0.0 
DO 15 K=l»4 
IF(K.NE.l)GO TO 1 
R1=R(I»J) 

Z1=Z(I»J) 

R2=R(I+1» J) 

Z2=Z(I+lr J) 

GO TO 4 

1 IF(K.NE.2)G0 TO 2 
R1=R(I+1»J) 

Z1=Z(I+1»J) 

R2=R(I+lr J+l) 

Z2=Z l 1+1 » J+l ) 

GO TO 4 

2 IF(K.NE.3)G0 TO 3 
R1=R(I+1»J+1) 

Z1=Z<I+1»J+1) 

R2=R<I»J+1) 

Z2=Z<I»J+1) 

GO TO 4 

3 IF (K«NE »4) GO TO 4 
R1=R C I » J+l ) 

Z1=ZU»J+1) 

R2=R(I» J) 

Z2=ZCI»J) 

4 R3=RC 
Z3=ZC 

AT=0.5*(R2*Z3-R3*Z2+R3*Z1-R1*Z3+R1*Z2-R2*Z1> 

A ( 1 » 1 ) = . 5* ( R2*Z3-R3*Z2 ) /AT 

A(2»1)=.5*(Z2-Z3)/AT 

A(3»1)=.5*(R3-R2)/AT 

A(1»2)=.5*(R3*Z1-R1*Z3)/AT 

A(2»2)=.5*(Z3-Z1)/AT 

A(3»2)=.5*(R1-R3)/AT 

A(1»3)=.5*(R1*Z2-R2*Z1)/AT 

A(2»3)=.5*(Z1-Z2)/AT 

A(3»3)=.5*(R2-R1)/AT 

DO 5 Ll=l»3 

GQ(L1)=0.0 

DO 5 L2=l»3 

DO 5 L3=l»2 

GT(Ll>L2>L3)=0.0 

5 CT(L1»L2»L3)=0.0 
IF(ABS(R1-R2) .LT.O.OOOOODGO TO 6 
AI_F12=(Z1-Z2)/(R1-R2> 

BET12=<Z2*R1-Z1*R2)/(R1-R2> 

GO TO 7 
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6 ALF12=0.0 
BET12=G.O 

7 IF(ABS(H2-R3) .LT. 0.000001)60 TO 8 
ALF23= ( Z2-Z3 ) / ( R2-R3 ) 

BET23= (Z3*R2-Z2*R3) / (R2-R3) 

GO TO 9 

8 ALF23=0. 

BET23=0. 

9 IF(ABS(R3-R1) .LT. 0.000001)60 TO 10 
ALF31=(Z3-Z1)/<R3-R1) 

BET31= <Z1*R3-Z3*R1 ) / ( R3-R1 ) 

60 TO li 

10 ALF31=0.0 
BET31=0.0 

11 CONTINUE 
RZ(1)=R1 
RZ(2)=R2 
RZ(3)=R3 
ZR(1)=Z1 
ZR(2)=Z2 
ZR(3)=Z3 

CALL INT6RL(RZ»ZR»X) 

GT(2»2'1)=X(3) 

6T(3»3»1)=X(3) 

GT(1»1»2)=X(1> 

6T(1»2»2)=X(2) 

GT(1»3»2)=X(6) 

6T(2»1»2)=X(2) 

GT(2»2»2)=X(3) 

GT(2»3»2)=X(7) 

6T(3»1»2)=X(6) 

6T(3»2»2)=X(7) 

6T(3*3»2)=X(10) 

DO 50 Ll=lr3 
DO 50 L2=l * 3 
DO 50 L3=l*2 

50 GT(L1»LZ»L3)=CK*GT(L1»L2»L3) 

IF(J.NE.JM1)G0 TO 51 
IF(K.NE.3)G0 TO 51 
SA=SQRT ( 1 . + ALF12**2 ) 

N=1 

GT ( 1 » 1 » N ) =6T (1'1'N) +CF*SA* < Rl**2-R2**2 ) /2 . 

GT(1»2»N)=GT(1»2#N) +CF*SA* (Rl**3-R2**3) /3. 0 
GT(l»3»N)=GT(l»3»N)+CF*SA*(ALF12*(Rl**3-R2**3)/3. 0+0.5 
1*BET12*(R1**2-R2**2) ) 

GT(2»1»N)=GT(1»2#N) 

GT(2>2'N)=GT(2>2rN)+CF*SA*0.25*(Rl**4-R2**4) 

GT(2»3»N)=GT(2»3»N)+CF*SA*(0.25*ALF12*(R1**4-R2**4) 

1 +BET12*(R1**3-R2**3)/3.Q) 

GT (3» 1 » N) =GT { 1 » 3 » N) 

GT(3»2»N)=GT(2»3»N) 

GT(3»3»N)=GT<3»3>N)+CF*SA*(0.25*ALF12**2*(R1**4-R2**4)+2.0*ALF12 
1 *BET 12* (Rl**3-R2**3) /3. 0+0 .5*BET12**2* (Rl**2-R2**2) ) 

GQ ( 1 ) =CF*SA*0 . 5* (Rl**2-R2**2) 

GQ(2)=CF*SA*(Rl**3-R2**3)/3.0 
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GQ ( 3) =CF*SA* ( ALF12* <R1**3-R2**3> /3. 0+0 .5*BET12* (Rl**2-R2**2) ) 
51 CONTINUE 

DO 48 Ll=l»3 

QR ( I ) =QR ( I ) +A < Ll » 2 ) *GQ ( Ll ) 

IF(I.EQ.IM)GO TO 48 

QR ( 1+1 ) =QR ( 1+1 ) +A ( Ll » 1 ) *GQ ( Ll ) 

48 CONTINUE 

DO 12 Ll=l t 3 
DO 12 L2=l»3 
DO 12 L3=l t 3 
DO 12 L4=l»3 
DO 12 L5=lr2 

12 CT(L1»L2»L5)=CT(L1»L2»L5)+A(L3»L1)*GT(L3»L4»L5)*A(L4»L2) 

DQ (K) =DQ (K ) +• 5*GT (2» 2» 2) *CP/CK 

IF(K.EQ.4)G0 TO 31 

DQ ( K+l ) =DQ ( K+l ) + . 5*GT ( 2 1 2 1 2 ) *CP/CK 

GO TO 32 

31 DQ(1)=DQ(1)+.5*GT(2»2»2)*CP/CK 

32 CONTINUE 

DO 14 Ll=l»2 

CQ(KrKrLl)=CQ(K#K»Ll)+CT(l»l»Ll) 

CQ(K»5+Ll)=CQ(Kr5»Ll)+CT(l»3»Ll) 

IF(K.EQ.4)G0 TO 30 
CQ(K»K+1*L1)=CT(1»2»L1) 

CQ(K+1»K+1»H)=CQ(K+1»K+1»L1)+CT(2»2»L1) 

CQ(K+1»5#L1)=CQ(K+1»5»L1)+CT(2»3»L1) 

GO TO 13 

oO CQ(1»1»L1)=CQ(1»1»L1)+CT(2»2»L1) 

CQ(1»4»H)=CT(1»2»L1) 

CQ(1»5»L1)=CQ(1»5»L1)+CT(2»3»L1) 

13 CQ(5»5»L1)=CQ(5»5»L1)+CT(3»3»L1) 

14 CONTINUE 
IF(J.NE.l)GO TO 15 
IF(K.NE.1)G0 TO 15 
IF ( I »EQ« IM) GO TO 15 
SAR=SQRT(1.+ALF12**2)/(R1-R2) 

QB ( I * I ) - SAR*(R2**2+Rl*R2-2.0*Rl**2)/6.0 
QB ( 1+1 » I ) = SAR* < 2 . 0*R2**2-R1*R2-R1**2 ) /6. 0 

15 CONTINUE 

DO 19 Ll=l»2 
DO 18 L2=2»5 
L4=L2-1 
DO 18 l3=1»L4 

18 CQ(L2»L3»Ll)=CQ(L3»L2»Ll) 

DO 19 L2=l»4 

DO 19 L3=l»4 

19 CQ(L2»L3>Ll)=CQ(L2»L3»Ll)-CQ(L2r5»Ll)*CQ(L3r5»Ll)/CQ(5»5»Ll) 

DO 40 L3=l»2 

Fl(J» J» l3)=F1( J» J»L3)+CQ(1»1»L3) 

Fl(J» J+1»L3)=F1 IJ» J+l»L3)+CQ(l»4rL3) 
Fl(J+l»J+l»L3)=Fl(J+l»J+l»L3)+CQ(4»4rL3) 

F2(vJ» J»L3)=F2( J> J»L3)+CQ(2»2»L3) 

F2(J» J+1»L3)=F2(J» J+1»L3)+CQ(2»3»L3) 

F2 ( J+l * J+l » L3 ) =F2 ( J+l > J+l » L3 ) +CQ ( 3 r 3 » L3 ) 

G( Jr J»L3)=G( J» J»L3)+CQ(1»2»L3) 

G(J» J+1»L3)=G(J» J+l»L3)+CQ(l»3rL3) 

G ( J+l » J» L3) =G (J+l » J» L3) +CQ ( 4» 2»L3) 
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G ( J+l t J4-1 , L3 ) =6 ( J+ 1 r J+l , L3 ) +CQ ( 4 1 3 1 l_3 ) 
40 CONTINUE 

Dl(J)=Dl(J)+DQU> 

D1(J+1)=D1( J+1)+DQ(4) 

D2(J)=D2(J)+DQ<2) 

D2(J+1)=D2(J+1>+DQ<3) 

21 CONTINUE 

DO 22 L3=l»2 
DO 22 L1=1»JM 
DO 22 L2=L1»JM 
IR=JM*(I-1) 

LR=IR+Ll 

LC=IR+L2 

C(LR»LC»L3)=F1(L1»L2»L3) 

22 F1(L1»L2#L3)=F2(L1»L2»U3) 

DO 26 L1=1»JM 
IR=JM*(I-1) 

LR=IR+Ll 
D (LR) =D1 (LI ) 

26 D1(L1)=D2(L1) 

IF(I.EQ.IM)GO TO 24 
DO 23 L3=l»2 

DO 23 L1=1»JM 
DO 23 L2=1»JM 
IR=JM*(I-1) 

LR=IR+U 

LC=IR+JM+L2 

23 C(LR»LC»L3)=G(L1»L2»L3) 

24 CONTINUE 

DO 27 L3=l » 2 
DO 27 Ll=2»M 
L1M=L1-1 
DO 27 L2=1»L1M 

27 C(Ll»L2»L3)=C(L2rLl#L3) 

RETURN 

END 
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ooo non odd 


STIFF 

SUBROUTINE STIFF (E» V»Al_F# IM> JM#N»CK»CF) 

COMMON K » Z 

COMMON/bTIF/KD»KS#PD»PS»PL»D 
REAL LAMfMU»KD»KT»KQ»KS 

DIMENSION R(15»5) >Z(15»5) »KD(15#15»15) >KS(15»15»15) »PD(15»5»15> » 
IPS ( 15» 5# 14) #PL(15»5»14) # D ( 75» 75) 

DIMENSION A ( 3» 3) t NN (10) r RT (9» 3 ) » KT (9 > 9) »KQ(15»15) » PT < 9» 3) » 
2PQ(15»5) #GT(3»3»2) »CQ<5»5*2) »CP(5»5) »DQ<4) »6Q(3) >CT(3»3r2) 
DIMENSION RZ(3) »ZR<3)*X(10) 

FN=N 

IM1=IM-1 

JM1=JM-1 

LAM=V*E/( (l.+V)*(l.-2.*V) ) 

MU=0.5*E/(1.+V) 

BET=ALF*(3.*LAM+2.*MU) 

El=LAM+2.*MU 

NN(1)=1 

NN(2)=3 

NN(3)=4 

NN(4)=6 

NN(5)=7 

NN(6)=9 

NN(7)=10 

NN(3)=12 

NN(9)=13 

NN(10)=15 


INITIALIZE 

DO 101 Ll=lrl5 
DO 101 L2=l»15 
DO 101 l3=1»15 
KS(L1»L2»L3)=0.0 

101 KD(L1»L2»L3)=0.0 
DO 102 Ll=l»l5 
DO 102 L2=l»5 
DO 102 L3=l»15 

102 PD(L1»L2»L3)=0.0 
DO 222 Ll=l»15 
DO 222 L2=l»5 
DO 222 L3=l»14 
PL (Ll »L2»L3) =0 • 0 

222 PS(Ll»L2rL3)=0.0 

OUTER LOOP ON I BEGINS HERE 

DO 402 I=1»IM1 

INNER LOOP ON J BEGINS HERE 

DO 21 J=1»JM1 
C INITIALIZE FOR I*J QUAD 
DO 104 Ll = l * 15 
DO 104 i_2=l»15 
104 KQ(LltL2)=0.0 
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DO 105 Ll=l t 15 
DO 105 L2=l»5 
105 PQ(L1»L2)=0.0 
R1=R(I»J) 

Z1=Z(I»J) 

R2=R(I+i» J) 

Z2=Z(I+1»J) 

R3=R(I+1» J+l) 

Z3=Z(I+1»J+1) 

R4=R(I» J+l) 

Z4=Z(I>J+1) 

R21=R2-R1 

R32=R3-R2 

R41=R4-R1 

R34=R3“R4 

Z21=Z2-Z1 

Z32=Z3-Z2 

Z41=Z4-Z1 

Z34=Z3-Z4 

AREA=.5*(R4l*Z41-R21*Z2l-R32*Z32+R34*Z34)-R32*Z21+R34*Z41 
AR=.5*(R41*Z4l)*(Rl+2.*R41/3.0)-.5*R21*Z2l*(Rl+2.*R2l/3. )-R32*Z2l* 
3(R2+.5*R32)-.5*R32*Z32*(R2+2.*R32/3. )+.5*R34*Z34*(R4+2.*R34/3. )+R3 
44*Z41*(R4+.5*R34) 

AZ;:.5*R41*Z41*<Zl+Z41/3. )-.5*R21*Z21*(Zl+Z21/3.)-R32*Z21*(Zl + .5*Z2 
51 ) 5*R32*Z32* (Z2+Z32/3. ) +. 5*R34*Z34* (Z4+Z34/3. ) +R34*Z41* (Zl+.5*Z4 
61 ) 

RC=AR/AREA 
ZC-AZ/ AREA 
DO 16 l_l=l,5 
DO 16 L2=l»5 
DO 16 L3=l»2 

16 CQ(L1»L2»L3)=0.0 
DO 17 Ll=l»4 

17 DQ (LI ) =0 • 0 
DO 15 K=l*4 
DO 601 Ll=l»9 
DO 601 L2=l»9 

601 KT(L1»L2)=0.0 
DO 602 Ll=l#9 
DO 602 L2=l»3 

602 PT(L1»L2)=0.0 
IF(K.NE.l)GO TO 1 
R1=R(I» J) 

Zl=Z(I»J) 

R2=R(I+1»J) 

Z2=Z(I+lr J) 

GO TO 4 

1 IF (K.NE.2) GO TO 2 
R1=R(I+1*J) 

Z1=Z(I+1» J) 

R2=R(I+1» J+l) 

Z2 = Z ( 1 + 1 » J+l ) 

GO TO 4 

2 IF(K.NE.3)G0 TO 3 
R1=R(I+1» J+l) 

Z1=Z(I+1» J+l) 

R2=R(I'J+1) 
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Z2=Z(If J+l) 

GO TO 4 

3 IF(K.NE.4)G0 TO 4 
Ri=R(I* J+l) 

Z1=Z<I»J+1> 

R2=R(I»J) 

Z2=Z<I»J) 

4 R3=RC 
Z3=ZC 

AT=0.5*IR2*Z3-R3*Z2+R3*Z1“R1*Z3+R1*Z2-R2*Z1) 

A ( 1 » 1 ) = . 5* ( R2*Z3-R3*Z2 ) /AT 

A(2»1)=.5*(Z2-Z3)/AT 

At3»l)=.5*(R3“R2)/AT 

A(1»2)=.5*(R3*Z1-R1*Z3)/AT 

A(2»2)=.5*(Z3-Z1)/AT 

A(3r2)=.5*(Rl“R3)/AT 

A ( 1 » 3 ) = . 5* (R1*Z2-R2*Z1 ) /AT 

A(2r3)=.5*CZl-Z2)/AT 

A(3»3)=.5*(R2-R1)/AT 

DO 5 Ll=l»3 

GQ(L1)=0.0 

DO 5 L2=l»3 

DO 5 L3=l#2 

GT(L1»L2»L3)=0.0 

5 CT(L1»L2#L3)=0.0 
IF(ABSCR1-R2) .LT. 0. 000001)60 TO 6 
ALF12=(Z1-Z2)/(R1-R2) 
BET12=(Z2*R1-Z1*R2)/(R1-R2) 

GO TO 7 

6 ALF12=0.0 
BET12=0.0 

7 IF(ABS(K2-R3) .LT. 0 . 000001 ) GO TO 8 
ALF23= ( Z2-Z3 ) / < R2-R3 ) 

BET23= (Z3*R2-Z2*R3) / (R2-R3) 

GO TO 9 

8 ALF23=0. 

BET23=0. 

9 IF(ABS(K3-R1) .LT.O.OOOOODGO TO 10 
ALF31=(Z3-Z1)/(R3-R1) 

BET31= CZ1*R3-Z3*R1 ) / ( R3-R1 ) 

GO TO 11 

10 ALF31=0.0 
BET31=0.0 

11 CONTINUE 
RZ(D=R1 
RZ(2)=R2 
RZ(3)=R3 
ZR(1)=Z1 
ZR(2)=Z2 
ZR(3)=Z3 

CALL INTGRL (RZ» ZR» X) 

GT(2»2»i)=X(3) 

GT(3»3»1)=X(3> 

GT(1»1»2)=X(1) 

GT(1»2»2)=X(2) 

GT(1»3»2)=X(6) 

GT(2»1»2)=X(2) 
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o o 


GT(2»2»2)=X(3) 

GT(2»3>«i)=X(7) 

GT(3»1»2)=X(6) 

GT(3»2»2)=X(7) 

GT (3»3»2)=X(10) 

R314=R3**4-R1**4 

R313=R3**3-R1**3 

R234=R2**4-R3**4 

R233=R2**3-R3**3 

R312=R3**2-R1**2 

R232=R2**2-R3**2 

G01=0 .25* ( ALF31-ALF12 ) *R314+ ( BET31-BET12 ) *R313/3. 0+0 • 25* ( ALF23- 
4ALF12) *R234+ (BET23-BET12) *R233/3. 0 
GD2=0 . 125* ( ALF31**2-ALFl2**2 ) *R314+ ( ALF31*BET31-ALF12*BET12) * R 

5313/3. 0+. 25* ( BET31**2“BET12**2 ) *R312+0. 125* (ALF23**2“ALF12**2)* 
6R234+(ALF23*BET23-ALF12*BET12)*R233/3.0+.25*(BET23**2-BET12**2)* 

2R232 

GD3=(ALF31**3-ALF12**3)*R314/12.0+(ALF31**2*BET31-ALF12**2*BET12 
6)*R313/3.0+0.5*(ALF31*BET31**2-ALF12*BET12**2)*R312+(BET31**3- BE 

7T12**3 ) * ( R3-R1 ) /3 • 0+ < ALF23**3-ALF12**3 ) *R234/ 12 • 0+ ( ALF23**2*BET23- 
SALF12**2*BET12 ) *R233/3 . 0+0 . 5* ( ALF23*BET23**2-ALF12*BET12**2 ) *R232+ 
9 (BET23** 3-BET12**3)*(R2“R3)/3.0 

FORM KT AND PT FOR TRIANGLE 
INSERT A 
Y 1=1 . 

Y2=l. 

IFCN.EQ.0)Y1=2. 

IF(N.EQ.0)Y2=0. 

DO 209 11=1 »3 
DO 209 Jl=l t 3 
AI=A<1»I1) 

BI=A(2»I1> 

DI=A(3»I1) 

AJ=A ( 1 » J1 ) 

BJ=AC2»J1) 

DJ=A(3»J1> 

GG1=GT(2»2»1) 

GG2=AJ*GT ( 1 » 2 » 2 ) +B J*GG1+D J*GT ( 2 1 3 r 2 ) 
GG3=AI*GT(l»2f2)+BI*GGl+DI*GT(2»3»2) 

GG4=AI *A J*GT ( 1 » 1 1 2 ) + ( A I*BJ+ A J*BI ) *GT ( 1 1 2 1 2 ) + < AI*DJ+A*J*DI ) *GT ( 
ll»3»2)+Bl*BJ*GGl+(Bl*DJ+BJ*Dl ) *GT(2»3»2)+DI*DJ*6T (3»3»2) 
GG5=AI*GG1+BI*GD1+DI*GD2 

GG6=A I * A J*GT ( 1 » 2 * 2 ) + ( A I *B J+A J*BI ) *GG1+ ( A I *DJ+A J*DI ) *GT ( 2 * 3 » 2 ) + 
1BI*BJ*GD1+ ( BI*DJ+BJ*DI ) *GD2+DI*D J*GD3 
GG7=AI*GG1+BI*GD1+DI*GD2 
DO 208 Ll=l * 3 
DO 208 L2=l » 3 
LR=(I1-1)*3+L1 
LC=('J1“1)*3+L2 
IF(Ll.NE.l) GO TO 202 
IF (L2.NE. 1 ) GO TO 200 
PT(LR* J1)=BET*Y1*GG5*BJ 

KT(LR»LC)=Y1*(E1*BI*BJ*GG1+LAM*(BI*GG2+BJ*GG3)+E1*GG4)+Y2*FN*FN 
2*GG4*MU+Y1*MU*DI*DJ*GG1 
GO TO 208 

200 IF (L2.NE.2) GO TO 201 
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KT ( LR f LC ) =Y 1 *FN* ( LAM*B I *GG2+E1*GG4 ) +Y2+MU* ( GG4-B J*GG3 ) *FN 
GO TO 208 

201 KT ( LR * LC ) =Y1* ( LAM*D J* ( BI*GG1+GG3) +MU*DI*B J*GG1 ) 

GO TO 208 

202 IF(L1.NE.2)G0 TO 205 
IF(L2.N£.1)G0 TO 203 
PT(LR» Jl)=-BET*FN*Y2*GG6 

KT ( LR f LC ) =Y1*FN* ( LAM*BJ*GG3+E1*GG4 ) +Y2*MU*FN* ( GG4-BI*GG2 ) 

GO TO 208 

203 IF(L2.NE.2)G0 TO 204 

KT(LRrLC)=Yl*FN*FN*El*GG4+Y2*MU*(Bl*BJ*SGl-Bl*GG2-BJ*GG3+GG4+Dl 

3*DJ*GG1) 

GO TO 208 

204 KT(LRrLC) =Y1*LAM*FN*DJ*GG3-Y2*MU*FN*DI *GG2 
GO TO 208 

205 IF (L2.NE • 1 ) GO TO 206 
PT(LRr Ji)=BET*Yl*DJ*GG7 

KT (LR» LC ) =Y1* (LAM* (DI*BJ*GG1+DI*GG2 ) +MU*DJ*BI*GG1 ) 

GO TO 208 

206 IF (L2.NE.2) GO TO 207 

KT ( LR » LC ) =Y1*LAM*FN*DI*GG2-Y2*MU*FN*DJ*GG3 
GO TO 208 

207 KT (LR t LC ) =Y1* ( E1*DI*DJ+MU*BI*BJ) *GG1+Y2*MU*FN*FN*GG4 

208 CONTINUE 

209 CONTINUE 

SUBROUTINE TO CALCULATE RT 

DO 210 I 1=1 f 9 
DO 210 Jl=l»3 

210 RT(I1*J1)=0.0 
Y3=l. 

IF(J.EQ.l.AND.K.EQ.l)GO TO 211 
IF(J.EQ.JM1.AND.K.EQ.3) GO TO 211 
GO TO 212 

211 CONTINUE 

GRl=R2**4/12.-R2**2*Rl**2/2.+2.*R2*Rl**3/3.-Rl**4/4. 

GR2=Rl**3*R2/6.-Rl**4/12.-Rl*R2**3/6.+R2**4/12. 

GR3=R2**4/4.-2.*Rl*R2**3/3t+Rl**2*R2**2/2.-Rl**4/12. 

BBB=BET*Y1*Y3/ (R2-R1 ) **3 

RT ( 1 » 1 ) =BBB* ( Z2-Z1 ) *GR1 

RT ( 3 » 1 ) =-l . *BBB* ( R2-R1 ) *GR1 

RT ( 4 » 1 ) =BBB* ( Z2-Z1 ) *GR2 

RT(6»1)=-1.*BBB*(R2-R1>*GR2 

RT(1»2)=RT(4»1) 

RT(3»2)=RT(6»D 

RT ( 4 > 2 ) =BBB* ( Z2-Z1 ) *GR3 

RT(6»2)=-1.*BBB*(R2“R1)*GR3 

DO 213 11=1 * 9 

DO 213 Jl=l»3 

213 PT(I1»J1)=PT(I1»J1)-RT(I1»J1) 

212 CONTINUE 

IF(I.EQ.1.AND.K.EQ.4) GO TO 221 
IF(I.EQ.IM1.AND.K.EQ.2)G0 TO 221 
GO TO 225 
221 CONTINUE 

GR4=Z2**4/12.-Z2**2*Zl**2/2.+2.*Z2*Zl**3/3.-Zl**4/4. 
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GR5=Zl**3*Z2/6.-Zl**4/12.-Zl*Z2**3/6.+Z2**4/12. 

GR6=Z2**4/4.-2.*Zl*Z2**3/3.+Zl**2*Z2**2/2.-Zl**4/12. 

GR7=Z2**3/3.-Z2**2*Zl+Z2*Zl**2-Zl**3/3. 

GR8=Zl**2*Z2/2.-Zl**3/6.-Zl*Z2**2/2.+Z2**3/6. 

GR9=Zl**2*Z2-Zl**3/3.-Zl*Z2**2+Z2**3/3. 

ALP=(R2-R1)/(Z2-Z1) 

BEP=(Rl*Z2-R2*Zl)/(Z2-Zl) 

BBB=BET*Y1*Y3/ (Z2-Z1 > **3 

RT ( 1 • 1 ) =BBB* ( Z2-Z1 ) * < ALP*GR4+BEP*GR7 ) 

RT (3* 1 ) =BBB* ( R2-R1 ) * < ALP*GR4+BEP*GR7 ) * ( -1 . 0 ) 

RT ( 4 » 1 ) =BBB* ( Z2-Z1 ) * < ALP*GR5+BEP*GR8 ) 

RT(6» 1 ) =BBB* ( R2-R1 ) * ( ALP*GR5+BEP*GR8 ) * ( - 1 . 0 ) 

RT(lr2)=RT(4fl) 

RT(3r2)=RT(6»l) 

RT ( 4 * 2 ) =BBB* ( Z2-Z1 ) * C ALP*GR6+BEP*GR9 ) 

RT ( 6 1 2 ) =BBB* ( R2-R1 ) * ( ALP*GR6+BEP*GR9 ) * ( -1 . 0 ) 

DO 223 Il=l»9 
DO 223 Jl=l»3 

223 PT ( 1 1 f Jl ) =PT ( 1 1 » J1 ) “RT ( 1 1 > J1 ) 

225 CONTINUE 

FORM CT 

DO 50 H=l»3 
DO 50 L2=l»3 
DO 50 L3=l»2 

50 GT(L1»L2»L3)=CK*GT(L1»L2»L3) 

IF ( J»NE. JM1 ) GO TO 51 
IF(K.NE.3)G0 TO 51 
SA=SQRT ( 1 . +ALF12**2 ) 

GT ( 1 » 1 » 1 ) =GT ( 1 » 1 t 1 ) +CF*SA* ( Rl**2-R2**2 ) /2 . 
GT(l»2»l)=GT<l»2»l)+CF*SA*(Rl**3-R2**3)/3. 

GT(l»3H)=GT(l»3»l)+CF*SA*(ALF12*(Rl**3-R2**3)/3.+.5*BET12*(Rl**2 
3”R2**2 ) ) 

GT(2»1»1)=GT(1>2»1) 

GT(2»2»1)=GT(2»2»1)+CF*SA* .25*<R1**4-R2**4) 

GT(2»3»i)=GT(2»3>l)+CF*SA*(0.25*ALF12*(Rl**4-R2**4)+BETl2*(Rl**3 
1 -R2**3)/3.0) 

GT(3»1»1) =GT ( 1 r 3 » 1 ) 

GT(3»2»1)=GT(2»3»1) 

GT(3»3»l)=GT(3>3»l)+CF*SA*<0.25*ALF12**2*(Rl**4-R2**4)+2.0*ALFl2 
1 *BETl2*(Rl**3-R2**3)/3.0+0.5*BET12**2*(Rl**2-R2**2) ) 

51 CONTINUE 

DO 12 Ll=l»3 
DO 12 L2=l»3 
DO 12 L3=l t 3 
DO 12 L4=l»3 
DO 12 L5=l»2 

12 CT(LlfL2»L5)=CT(LlfL2»L5)+A(L3»Ll)*GT(L3rL4rL5)*A(L4fL2) 

NOW FOR THE QUADRILATERAL 

DO 300 Kl=lr3 
DO 300 K2=l>3 
KR=3*(K-1)+K1 
KC=3*(K-1)+K2 
KR5=12+K1 

- 78 - 



KC5=12+K2 

KRR=3*K+K1 

KCC=3*K+K2 

KQ(KR»KC)=KQ(KR»KC)+KT(K1»K2) 

KQ<KR»KC5)=KQ(KR»KG5)+KT(Kl»K2+6) 

KQ(KR5>KC)=KQ(KR5>KC)+KT<K1+6»K2) 

IF(K.EQ.4) GO TO 301 

KQ(KR»KCC)=KQ(KR»KCC)+KT(Kl»K2+3) 

KQ(KRR»KC)=KQ(KRR»KC)+KT(K1+3#K2) 

KQ(KRRrKCC)=KQ(KRR»KCC)+KT(Kl+3»K2+3) 

KQ(KRR»KC5)=KQ(KRR»KC5)+KT(Kl+3»K2+6) 

KQ(KR5»KCC)=KQ(KR5»KCC)+KT(Kl+6#K2+3) 

GO TO 302 

301 KQ(K1*K2)=KQ(K1»K2)+KT (Kl + 3 * K2+3 ) 
KQ(Kl»K2+9)=KQ(Kl»K2+9)+KT(Kl+3»K2) 
KQ(Kl+9»K2)=KQ(Kl+9*K2)+KT<Kl»K2+3) 
KQ(KRR»K2)=KQ(KRR»K2)+KT<Kl+6»K2+3) 
KQ(Kl»KCC)=KQ(KlrKCC)+KT(Kl+3»K2+6) 

302 KQ(KR5»KC5)=KQ(KR5»KC5)+KT(Kl+6»K2+6) 

300 CONTINUE 

DO 303 Kl=lr3 
KR=3*(K-1)+K1 

PQ(KR»K)=PQ(KRrK)+PT(Kl,l) 
PQ(KR»5J=PQ(KR»5)+PT(K1»3) 
PQ(K1+12»K)=PQ(K1+12»K)+PT(K1+6»1) 
IF(K.EQ.4) GO TO 304 
PQ(KR»K+1)=PQ(KR»K+1)+PT(K1»2) 

PQ (KR+3»K) =PQ (KR+3» K) +PT (Kl+3» 1 ) 

PQ ( KR+3 » K+l ) =PQ ( KR+3 » K+l ) +PT < Kl+3> 2 > 

PQ (KR+3 » 5) =PQ (KR+3 » 5) +PT (Kl+3» 3) 
PQ(K1+12»K+1)=PQ(K1+12»K+1)+PT(K1+6»2) 

GO TO 305 

304 PQ(K1»4)=PQ(K1»4)+PT(K1+3»1> 
PQ(Kl»5)=PQ(Kl»5)+PT(Kl+3»3) 
PQ(K1+9»1)=PQ(K1+9»1)+PT(K1»2) 

PQ (K1 » 1 ) =PQ (Kl » 1 ) +PT (Kl+3» 2 ) 

PQ(Kl+l2»l) =PQ (Kl+12» 1 ) +PT ( Kl+6» 2) 

305 PQ(Kl+l2»5)=PQ(Kl+12»5)+PT(Kl+6»3) 

303 CONTINUE 

DO 14 Ll=l»2 

CQ(K»K»Ll)=CQ(K»K»Ll)+CT(lrl»Ll> 
CQ(K»5»L1)=CQ(K»5»L1)+CTC1»3»1-1) 
IF(K.EQ.4)G0 TO 30 
CQ(K»K+1*L1)=CT(1»2*L1) 

CQ(K+l»K+lfLl)=CQ(K+l*K+l»Ll)+CT(2»2»Ll> 

CQ(K+l»5»Ll)=CQ(K+l»5»Ll)+CT(2r3»Ll) 



GO TO 13 




2»L1 ) 

30 

CQ ( 1 ’ 1 ’ Ll ) 

=CQ( 

1 r 1 » Ll ) +CT 

(2 » 


CQ(1»4»L1) 

=CT( 

l»2rLl) 




CQ(1»5»L1) 

=CQ( 

1 » 5» Ll ) +CT 

(2» 

3»L1) 

13 

CQ(5»5»Ul) 

=CQ( 

5»5»L1)+CT 

(3» 

3»L1) 

14 

CONTINUE 






DO 18 Ll=l 

1 2 





DO 19 L2=2 

1 5 





L4=L2-1 






DO 19 L3=l 

»L4 





19 CQ(L2>L3»L1)=CQ(L3»L2»L1) 
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18 CONTINUE 
15 CONTINUE 

DO 500 Ll=l»5 
DO 500 L2=l»5 

500 CP(Ll»L2)=C<si(Ll »L2»1) +FN*FN*CQ < LI » L2 » 2 > 

C 

c need to eliminate middle node 

c 

C INSERT A 

Nl=3 

IF(N.NE.O)GO TO 801 
Nl=2 

DO 800 Ml=l t 10 
M3=NN(M1) 

DO 802 M5=l»5 
802 PQ(M1»M5)=PQ(M3»M5) 

DO 800 M2=l»10 
M4=NN(M2) 

800 KQ(M1»M2)=KQ(M3»M4) 

801 CONTINUE 
N2=2*N1 
N3=3*Nl 
N4=4*Nl 

DO 306 K1=1»N1 
DO 306 K2=1#N1 

306 D ( K1 » K2 ) =KQ ( K1+N4 1 K2+N4 ) 

CALL INVERT (DrNl»75) 

DO 307 K1=1»N4 
DO 307 K2=lrN4 
DO 307 K3=lrNl 
DO 307 K4=1»N1 
L3=K3+N4 
L4=K4+N4 

307 KQ(Kl»K2)=KQ(Kl»K2)-KQ(Kl»L3)*D(K3rK4)*KQ(L4»K2) 

DO 309 Ll=l t N4 
DO 309 L2=l»4 

309 PQ(Ll»L2)=P0(Ll»L2)-PQ<Ll»5)*CP(5rL2)/CP(5»5) 

DO 308 L1=1»N1 

DO 308 L2=l»4 

308 PQ(N4+Ll»L2)=PQ(N4+Ll»L2)“PQ<N4+Llr5)*CP(5»L2)/CP(5r5) 
DO 310 L1=1»N4 

DO 310 L2=l»4 
DO 310 K1=1»N1 
DO 310 K2=1»N1 

310 PQ(Ll»L2)=PQ(Ll»L2)-KQ(Ll»Kl+N4)*D(KlrK2)*PQ(N4+K2»L2) 
C 

C ASSEMBLE THE ROW MATRICES KD»KS»PD»PS 
C 

DO 400 K1=1»N1 
DO 400 K2 = l * Nl 
KR=N1*(J-1)+K1 
KC=N1*(U-1)+K2 
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KD(KRrKCr I)=KD(KRrKCrI)+KQ(KlrK2> 

KD(KRrKC+Nlr I)=KD(KRrKC+Nlr I)+KQ(KlrK2+N3) 

KD (KR+Nl rKCr I )=KD( KR+Nl rKCr I ) +KQ(Kl+N3rK2) 

KO (KR+Nl r KC+N1 r I ) =KD ( KR+Nl r KC+Nl r I ) +KQ ( K1+N3 » K2+N3 ) 

KD (KRr KC rI+l)=KD(KRrKCrI+l)+KQ (Kl+Nl rK2+Nl> 

KO ( KR r KC+N1 r 1+1 ) =KD ( KR r KC+Nl r 1+1 ) +KQ ( Kl+Nl * K2+N2 ) 

KD(KR+NlrKCr I+1>=KD (KR+Nl »KCr I+l)+KQ(Kl+N2rK2+Nl) 

KD (KR+Ni r KC+Nl r 1+1 ) =KD ( KR+Nl r KC+Nl r 1+1 ) +KQ ( K1+N2 r K2+N2 ) 
KS(KRrKCr I)=KS(KRrKCrI)+KQ(KlrK2+Nl) 

KS (KRr KC+Nl r I ) =KS (KRr KC+Nl rI)+KQ(Klr K2+N2) 

KS (KR+Nl r KC r I ) =KS (KR+Nl »KC r 1 ) +KQ (Kl+N3r K2+N1 ) 

KS (KR+Nl r KC+Nl * I > =KS (KR+Nl r KC+Nl r 1 ) +KQ ( K1+N3 r K2+N2 ) 

400 CONTINUE 

DO 401 Kl=lrNl 
KR=N1*(J-1)+K1 

PD(KRr Jr I)=PD(KRr Jr I ) +PQ (K1 r 1 ) 

PD(KR» J+lrI)=PD(KRr J+lrI)+PQ(Klr4) 

PD(KR+Nlr Jr I)=PD(KR+Nlr Jr I)+PQ(K1+N3rl) 

PD (KR+Nl r J+lr I)=PD(KR+N1» J+lr I ) +PQ (Kl+N3r 4) 

PD (KRr Jr 1+1) =PD (KRr Jr I + D +PQ( Kl+Nl r 2) 

PD (KR r J+l r 1 + 1 ) =PD (KRr J+l r 1 + 1 ) +PQ( Kl+Nl r 3) 

PD (KR+Nl » Jr I+1)=PD (KR+Nl #J» I+1)+PQ(K1+N2r2) 

PD ( KR+Nl r J+lr 1+1 ) =PD (KR+Nl r J+lr I+1)+PQ(K1+N2r3) 

PS (KR r J r I ) —PS (KRrJrI)+PQ(Kl»2) 

PS(KRr J+lrI>=PS(KRr J+l r I ) +PQ ( K1 r 3) 

PS ( KR+Nl r Jr I )=PS( KR+Nl r J r I ) +PQ ( K 1+N3 r 2 ) 

PS ( KR+Nl r J+l r I ) =PS ( KR+Nl r J+l r I ) +PQ ( K1+N3 r 3 ) 

PL(KRr JrI)=PL(KRr J r I ) +PQ ( Kl+Nl r 1 ) 

PL (KRr J+l r I ) =PL(KRr J+l r I ) +PQ ( Kl+Nl r 4 ) 

PL ( KR+Nl r J r I ) =PL ( KR+Nl r J r I ) +PQ ( K1+N2 r 1 ) 

PL ( KR+Nl r J+l r I ) =PL ( KR+Nl r J+l r I ) +PQ ( K1+N2 r 4 ) 

401 CONTINUE 
21 CONTINUE 

402 CONTINUE 
RETURN 
END 
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INTGRL 

SUBROUTINE INTGRL (R » Z» X ) 

REAL IC0N»IC0NP#IZ»IZP#IZ2#IZ2P 
DIMENSION R(3) » Z ( 3 ) » X ( 10 > > XI < 10 > »AI(10) 

RI=R(1> 

RJ=R(2) 

RK=R(3) 

DATA (XI ( I ) # AI ( I ) » 1=1 » 10 ) /-. 97390653* . 066671344# 86506337# . 1494513 
l5»-.67940957» .21908636#-. 43339539# .26926672 #-.14887434# .29552422# 
2.14887434# .29552422# .43339539# .26926672# .67940957# .21908636, 
3.86506337# .14945135# .97390653# .066671344/ 

DO 2001 Nl=l#10 
2001 X(N1)=0. 

C**** CALCULATION OF INTEGRALS BY GAUSSIAN QUADRATURE 
70 RMIN=AMIN1(RI»RU»RK) 

RMAX=AMAX1(RI#RU#RK) 

DO 7 Nl=l#3 

7 IF ( ABS (R (N1 ) -RMIN) .LE. 0 . 00001 ) I1=N1 
DO 8 Nl=l # 3 

8 IF ( ABS (R (N1 ) -RMAX ) .LE. 0.00001) I3=N1 
DO 9 Nl=l»3 

9 IF(N1.NE.I1.AND.N1.NE.I3) I2=N1 

R1=R(I1) 

R2=R(I2) 

R3=R(I3) 

Z1=Z(I1) 

Z2=Z(I2) 

Z3=Z( 13) 

FAC=1.0 

DR12=ABS(R1-R2) 

DR13=ABS(R1-R3) 

IF (Rl.GT. 0.0001) GO TO 100 

IF (DR12.LT. 0.0001. OR. DR13.LT.0. 0001 )FAC=1000.0 
100 CONTINUE 

S12=(Z2-Z1)/(R2-R1) 

S13=(Z3-Z1)/(R3-R1) 

S23 = ( Z3-Z2 ) / ( R3-R2 > 

DR=R2-R1 

DRP=R3-R2 

DO 12 Nl=l»10 

RR=Rl+DR*(XI(Nl)+l.)/2. 

RRP=R2+0RP* ( X I ( N1 ) +1 . ) /2 . 

ZZ1=S13*(RR-R1)+Z1 
ZZ1P=S13* ( RRP-R1 ) +Z1 
ZZ2=S12*(RR-Rl)+Zl 
ZZ3=S23* ( RRP-R2 ) +Z2 
IC0N=ABS(ZZ2-ZZ1) 

IC0NP=ABS(ZZ3-ZZ1P) 

lZ=(ZZl**2-ZZ2**2)/2. 

IF (ZZl .LT • ZZ2) IZ= -IZ 
IZP=<ZZlP**2-ZZ3**2)/2. 

IF (ZZlP.LT .ZZ3) IZP= -IZP 
IZ2=ABS (ZZ2**3-ZZ1**3) /3. 
lZ2P=ABb(ZZ3**3-ZZlP**3)/3. 

DO 10 N2=l#5 

X ( N2 ) =X ( N2 ) +A I ( N1 ) * I CONP*RRP** ( N2-2 ) *DRP 

10 IF(ABS(RR) .GT. 0.0000001) 
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1 X ( N2 ) =X l N2 ) + A I < N1 ) * ( ICON*RR** < N2-2 ) *DR ) 

DO 11 N*=6»9 

X (N2) =X (N2) +AI (N1 ) *IZP*RRP** ( N2-7) *DRP 

11 IF(ABS(RR) .GT. 0.0000001) 

IX ( N2 ) =X ( N2 ) + A I < N1 ) * ( I Z*RR** ( N2-7 ) *DR ) 

12 X(10)=X(10)+AI <N1)*(IZ2/RR*DR+IZ2P/RRP*DRP) 
DO 13 Nl=l t 10 

13 X(Nl)=X(Nl)/2. 

X ( 1 ) =FAC*X ( 1 ) 

X(6)=FAC*X(6) 

X(10)=FAC*X(10) 

RETURN 

END 
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SUBROUTINE POTTER ( A* B» C r DL» Z» MAT* LCV * AS *DK* DKK* DKKK* IPIVOT* INDEX* 
1ZZ»SPP»SPPP»Q»P*SAP»SS*SP*SSS»MM»LL) 

CHANGED 5 MARCH 71 FOR DIMENSION STATEMENT CONTINUITY 
C* VARIABLES MM AND LL ADDED TO ARGUMENT LIST 

DIMENSION AS(LL*LL) »DK(LL*LL) »DKK(LL*LL) »DKKK (LL* LL) » 

1 ZZ(LL*1) »SPP(LL) »SPPP(LL) »Q(MM»LL) *Z(MM»LL) *P(MM»LL»LL) * 

2 SAP ( LL ) » SS ( LL ) * SP ( LL ♦ LL ) » SSS C LL ) * A ( MM » LL * LL ) * B ( MM * LL » LL ) * 

3 C ( MM * LL » LL ) * DL ( MM » LL ) 

DOUBLE PRECISION DETERM 
DIMENSION IPIVOT (LL) » INDEX (LL* 2) 

M=MAT 


C 

C 

C 


N=M-1 
NMAX=LL 
DO 1 I=1*LCV 
DO 1 J=1*LCV 
1 C(M*I*U)=0. 

LOGIC TO STATEMENT 22 CALCULATES PI AND Q2 MATRICES 


DO 4 I=1*LCV 
DO 4 J=1*LCV 
4 ASd* J)=A(2*I* J) 

CALL MATINV (AS* LCV »ZZ*0»DETERM» IPIVOT* INDEX 


NMAX* ISCALE) 


DO 15 1=1 * LCV 
DO 15 J=1*LCV 
DK(I*J)=.0 
DO 15 K=1 » LCV 

15 DK ( I » J ) =DK (I*J)+B(1»I*K) *AS ( K » J ) 

DO 16 1=1 *LCV 

DO 16 J=i * LCV 
DKK ( I * J) =0 • 

DO 16 K=1 » LCV 

16 DKK ( I * U ) =DK ( I * K ) *B ( 2 * K * U ) +DKK(I*J) 

DO 17 1=1 *LCV 

DO 17 J=1 *LCV 

17 DKKK ( I » J) =DKK ( I » J) “C ( 1 » I * J) 

CALL MATINV (DKKK » LCV > ZZ» 0 » DETERM* IPIVOT * INDEX*NMAX» ISCALE) 

DO 18 1=1 * LCV 
DO 18 J=1 »LCV 
SP(I*J)=0. 

DO 18 K=1 * LCV 

18 SP ( I » J > =SP ( I * J) +DK (I*K)*C(2*K»U) 

DO 19 I=1»LCV 

DO 19 J=1»LCV 
P(2»I*J)=0. 

DO 19 K=1 »LCV 

19 P(2»I»J)=P(2*I*J)+DKKK(I*K)*SP(K*J) 

DO 20 1=1 *LCV 

SPP(I)=G. 

DO 20 K=1 »LCV 

20 SPP ( I ) = SPP ( I ) +DK ( I * K ) *DL ( 2 * K ) 

DO 21 1=1 » LCV 

21 SPPP(I)=SPP(I)-DL(1*I) 

DO 22 1=1 »LCV 
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OOO OOO .0.000 


Q(2*I)=G. 

DO 22 K=1»LCV 

22 Q(2»I)=Q(2»1)+DKKK(I*K)*SPPP<K) 

LOGIC TO STATEMENT 23 CALCULATES* IN A LOOP* THE Q3....QM AND THE 
P2 TO PM MATRICES 

DO 23 L=3»M 
DO 24 I=1*LCV 
DO 24 J=1*LCV 
DK(I*J)=0. 

DO 25 K=1 » LCV 

25 DK ( I * J ) =DK ( I * J ) +A ( L * I » K ) *P ( L-l » K » J ) 

24 DK(I*J)=-DK(I*J)+B(L»I*J> 

CALL MATINV (DK* LCV * ZZ * 0 #DETERM» IPI VOT» INDEX»NMAX* ISCALE) 

DO 26 I=1*LCV 
DO 26 J=1»LCV 
P(L»I*J)=0. 

DO 26 K=1*LCV 

26 P(L»I*J)=P(L»I*J)+DK(I»K)*C(L»K»J) 

DO 27 1=1 »LCV 

SAP(I)=0. 

DO 28 K=1 *LCV 

28 SAP(I)=SAP(I)+A(L»I»K)*Q(L-1»K) 

27 SaP(I)=-SAP(I)+DL(L»I) 

DO 29 1=1 » LCV 
Q(L»I)=0. 

DO 29 J=1 » LCV 

29 Q(L»I)=0(L»1)+DK(I»J)*SaP( J) 

IF(L.EQ.M) GO TO 50 

GO TO 23 

50 DO 30 1=1 *LCV 

30 Z(L»I)=Q(L*I) 

23 CONTINUE 

LOGIC TO STATEMENT 32 CALCULATES* IN A LOOPr THE Z3 TO Z(M-l) MATRICES 
M=N-2 

DO 32 L=1*M 
K=M-L+3 
DO 32 1=1 » LCV 
Z(K»I)=0. 

DO 31 J=1 » LCV 

31 Z(K*I)=Z(K»I)+P(K»I*J)*Z(K+1*J) 

32 Z (K» I ) =~Z (K» I ) +0 (K * I ) 

REMAINING LOGIC CALCULATES zi AND Z2 

DO 33 1=1 » LCV 
SAP ( I ) =• 0 
DO 33 J=1 » LCV 

33 SAP(I)=SAP(I)-SP(I» J)*Z(3* J) 

DO 34 1 = 1 * LCV 

34 SAP C X > =SAP ( I ) +SPPP ( I ) 

DO 35 1 = 1 » LCV 
Z(2»I)=.0 
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DO 35 J=1 » LCV 

o5 Z(2»I)=Z(2»I) +DKKK ( I » J) *SAP ( J) 
DO 36 1=1 > LCV 
SS(I)=.0 
DO 36 J=1 »LCV 

36 SSII)=SS(I)+B(2»I»J)*Z(2»J) 

DO 37 1=1 rLCV 

SSS(I)=.0 
DO 37 J=1»LCV 

37 SSS(I)=SSS(I)+C(2> It J)*Z(3» J) 
DO 38 1=1 » LCV 

38 SSSCI)=-SSS(I)-SS(I)+DL(2rI) 

DO 40 1=1 »LCV 

Z(1»I)=.0 
DO 40 J=1 #LCV 

40 Zll»I)=Z(l»I) +AS ( I » J) *SSS ( J) 
RETURN 
END 


- 86 - 



INVERT 

SUBROUTINE INVERT (D» ACT #DIM) 

C INVERSION OF SYMMETRIC MATRIX 
INTEGER ACT»DIM 
DIMENSION D(DIM»DIM) »L0C(76> 

DOUBLE PRECISION DP 

DP=1.D0 

DO 1 N=1 t ACT 

1 LOC(N)=N 

DO 6 Nl=l t ACT 
M=0 

PIVOT=0. 

DO 2 N2=Nlr ACT 
NN=L0C<N2) 

IF (ABS(D(NN»NN) ) .LE.ABS(PIVOT) ) GO TO 2 
M=N2 

PIVOTrD (NNrNN) 

2 CONTINUE 

IF (M.EQ.O) GO TO 8 
N=LOC(M) 

LOC(M)=LOC(Nl) 

LOC(Nl)=N 

D(N»N)=-1. 

DO 3 J=1 » ACT 

3 D(Nr J)=D(N»J)/PIVOT 
DO 5 I1=1»ACT 
I=LOC(Il) 

IF (N.EQ.I.OR.D(I»N) .EQ.O.) GO TO 5 

DO 4 J1=I1,ACT 

J=LOC(Jl) 

IF (N.EQ.J) GO TO 4 

D(I» J)=U(I» J)-D(I»N)*D(n»U)*DP 

D ( J» I ) =0 ( I » J ) 

4 CONTINUE 

5 CONTINUE 

DO 6 1 = 1 » ACT 

6 D(I»N)=U(N»I) 

DO 7 1 = 1 » ACT 
DO 7 J=i » ACT 

7 D(IrJ)=-D(I»J) 

RETURN 

8 WRITE(6»9) 

9 FORMAT (42H0MATRIX IS SINGULAR - EXECUTION TERMINATED ) 
STOP 

END 
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c 

c 

c 


c 

c 

c 


MATINV 

SUBROUTINE MATINV (A»N»B»M»DETERM» IPIVOTr INDEX »NMAX’ I SC ALE) 

matrix inversion with accompanying solution of linear equations 

DIMENSION AlNMAX>N)»B<NMAX»M)»IPIV0T(N)’INDEX(NMAX»2) 
EQUIVALENCE (IROW’JROW), ( ICOLUM’ JCOLUM) » ( AMAX» T ’ SWAP) 

INITIALIZATION 

b ISCALE=U 
b Rl=10. 0**18 
7 R2=1.0/K1 
10 DETERM=1.0 
lb DO 20 J=1»N 
20 IPIVOT(J)=0 
30 DO 550 1=1 » N 


C 

C 

C 


c 

c 

c 


SEARCH FOR PIVOT ELEMENT 


40 

4b 

50 

60 

70 

BO 

8b 

90 

9b 

100 

10b 

10b 


110 


IF 

DO 

IF 

IF 


60 » 105’ 60 


AMAX = 0 • 0 
DO 105 J=1»N 

(IPIVOT(J)-l) 

100 K=1’N 

(IPlVOT(K)-l) 80 » 100’ 740 
(ABS(AMAX)-ABS(A(J’K) ) )85’100»100 
IR0W=J 
1C0LUM=K 
AMAX=A( J»K) 

CONTINUE 
CONTiNUt 

IF (AMAx) H0’106»110 
DETERM=0.0 

iscale=o 

GO TO 740 

IPIVOT ( ICOLUM ) =IPIV0T ( ICOLUM) +1 

interchange rows to put pivot element on diagonal 

130 IF (IROW-ICOLUM) 140’ 2b0» 140 
140 DETERM=“DETERM 
150 DO 200 L=1»N 
1 60 SWAP=A( IROW’L) 

170 A ( IROW ’ i_) = A ( ICOLUM ’L) 

200 A(ICOLUM’L)=SWAP 
20b IF (M) 260’ 260’ 210 
210 DO 250 L=l» M 
220 SWAP=B< IROW’L) 

230 B(IROW»D=B(ICOLUM»L) 

250 B(ICOLUM’L)=SWAP 
260 INDEX(I»1)=IR0W 
270 INDEX ( I »2) = ICOLUM 
310 PIVOT=A( ICOLUM’ ICOLUM) 


C 

C 

C 


SCALE THE DETERMINANT 
xOOO PIVOTI=PIVOT 
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1005 IF ( ABS ( DETERM) -Rl> 1030 #1010 #1010 
1010 DETERM=UETERM/R1 
lSCAL£=iSCALE+l 

IF ( ABS ( DETERM) -R1 > 1060 » 1020 # 1020 
l020 DETERM=DETERM/R1 
ISCALE=1SCALE+1 
BO TO 1060 

1030 IF ( ABS (UETERM >-R2) 1040 » 1040 » 1060 
1040 DETERM=UETERM*Ri 
1SCALE=1SCALE-1 

IF ( ABS ( UETERM )-R2) 1050 #1050 #1060 
1050 DETERM=UETERM*R1 

ISCALE=ISCALE-1 . 

1060 IF ( ABS (PIVOTI)-Rl) 1090 #1070 #1070 
1070 PIV0TI=PIV0TI/R1 
ISCALE=ISCALE+1 

1F(ABS(PIVOTI)-Rl)320#l0ti0»l080 

1080 PIV0TI=PIV0TI/R1 
1 1SCALE=ISCALE+1 
GO TO 320 

1090 IF (ABS (PIVOT I ) -R2) 2000 # 2000 » 320 
*000 PIV0TI=PIV0TI*R1 
ISCALE=ISCALE-1 

IF (ABS (PIVOT I )“R2) 2010 #2010 #320 
2010 PIVOTI=PIVOTI*Rl 
ISCALE=1SCALE-1 
320 D£TERM=UETERM*PIVOTI 


DIVIDE PIVOT ROW BY PIVOT ELEMENT 

330 A ( ICOLUM# ICOLUM) =1 • 0 
340 DO 350 l=1#N 

350 A ( ICOLUM »L)=A( ICOLUM# L) /PIVOT 
355 IF ( M) 3ti0» 300# 360 
360 DO 370 L=1#M 

370 B ( ICOLUM# L ) =B ( ICOLUM# L) /PIVOT 


C 

C 

C 


REDUCE IMON-PIVOT ROWS 


300 

390 

400 

420 

430 

450 

45b 


DO 550 L1=1»N 

IF (Ll-ICOLUM) 400# 550# 400 
T=A(L1» ICOLUM) 

A (LI # ICOLUM ) =0 • 0 
DO 450 L=1#N 

A(Ll#L)=A(Li»L)-A(ICOLUM»L)*T 

1F(M> 550# 550# 460 


460 DO 500 L=1»M 

500 B (LI »L)=B (LI »L)“B( ICOLUM# L)*T 
550 CONTINUE 


INTERCHANGE COLUMNS 

600 DO 710 I=1#N 
610 L=N+1-I 

620 IF (INDEX(L»1)-INDEX<L»2>> 
630 JR0W=INUEX(L»1> 

640 JC0LUM=INDEX(L»2) 


630# 710# 630 
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650 00 705 K=1»N 
660 SWAP=A<K» JRQW) 

670 A(K» JR0W)=A<K» JCOLUM) 
700 A (K* JCOlUM) =SWAP 
70b CONTINUE 
710 CONTINUE 
740 RETURN 

end 



MAIN 

C PROGRAM TO COMPUTE OPTIMAL CONTROL 

C 

c 

DIMENSION I HP ( 50 ) » ISP ( 180 ) # A ( 180 » 169) t B < 50 » 50 ) » C ( 50 r 180 ) » D( 50 » 50 ) 
1 WD 1 1QU ) » WW ( 180 ) 

bn = i • o 

C INITIALIZATION 

DO 12 1 = 1 » 180 
WD ( I ) =0 • 0 
WW(I)=0.0 
ISP(I)=U 
DO 11 J=i t 50 

11 C(J»I)=0.0 

DO 12 J=1 t 168 

12 A(I»J)=0.0 
DO 20 1=1 r 50 
IHP( I ) =0 • 0 
DO 20 J=1 r 50 

20 B(I»J)=0.0 

READ l 5» 50) !MrKM»DI»DO»FNO 
50 FORMAT (2I5»3F10«5) 

RI=Dl/2.0 
R0=D0/2 • 0 
C 

C INPUT SAMPLE POINT AND HEATER LOCATIONS 

C 

READ (5» 100 ) NHP»NSP 
100 FORMAT (215) 

READ (5» 110 ) (IHP(I) »I=1»NHP) 

110 FORMAT ( 1615) 

READ (5» 120 ) (ISP(I) » I=1»NSP) 

120 FORMAT ( 1615) 

WRITE(6»259) 

259 FORMAT (29H INDICES OF HEAT PATCH P0lNTS//4Xr 1HI r2X»3HlHP) 

DO 260 L1=1»NHP 

260 WRITE(6.261) L1»IHP(L1) 

261 F0RMAT(2I5) 

WRlTE(6r 262) 

262 FORMAT (25H INDICES OF SAMPLE POINTS//4Xr 1HI »2Xf 3HISP) 

DO 263 Ll=l * NSP 

263 WRITE(6»261) L1»ISP(L1) 

C 

C INPUT MATRIX A AND COMPUTE REDUCED MATRIX 

C 

READ (4 * 130) ( (A(I* J) » 1=1 » 180) »J=1»169) 

130 FORMAT (6E13. 8) 

C 

C SUBTRACT OUT SPHERICAL PART OF A(I»J) 

C 

DIMENSION THETU5) #BE<15) »Y<15> 

FIM1=IM-1 

fkm=km 

RF=2»0*DO*FNO 
SI=.5*DI/RF 
C I=SQRT ( 1 « 0”SI*SI ) 

SO=.5*DO/RF 
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non 


CO=SQRT(1.0-SO*SO) 
THETI=ATAN(5I/CI) 
THETO=ATAN(SO/CO) 

DTHET= ( THETO-THETI ) /FlMl 

theto=theto-dthet 

B2-0 * 0 

DO 39 1=1 # IM 

THET(I)=THETI+(1-1)*DTHET 
39 BE(I)=COS(THET(I) ) /COS ( THETO) -1 . 0 
DO 63 KH=1»NHP 
JH=IHP(KH) 

WB=0.0 
B2=0.0 

DO 61 16=1 t NSP 
I=ISP(IS)/IM 
I=ISP(IS)-I*IM 
JS=ISP(1S) 

WB=WB+A l JS ' JH ) *BE ( I ) 

61 B2=B2+BE(I)**2 
DO 62 1=1 » IM 

62 Y(I)=WB/B2*BE(I) 

DO 63 IS=1 r NSP 
I=ISP(IS)/IM 
I=ISP(IS)-I*IM 
JS=ISPCIS) 

63 A(JS»JH)=A(JS»JH)-Y(I) 

DO 150 I=1»NSP 
DO 150 J=1»NHP 
IS=ISP(I) 

JH=IHP(J) 

150 A ( I * J) =A ( IS» JH) 

TO COMPUTE I-A(1/ATA)aT 

DO 200 i=l»NHP 
DO 200 J=lrNHP 
DO 200 K=1»NSP 

200 B(I» J)=B(I» J)+A(K»I)*A(K» J)*BN 
DO 198 L1=1»NHP 
DO 198 L2=1»NHP 

198 D(Ll»L2)=B(LlrL2)/BN 
call invert (b»nhp»50) 

DO 199 L1=1»NHP 
DO 199 L2=1»NHP 

199 B(L1»L2)= BN *B(L1»L2) 

DO 197 I=1#NHP 
DO 197 J=1»NHP 
DO 197 K=1»NHP 

197 C(I»J)=C(I» J)+B(I»K)*D(K»J) 

DO 195 I=1»NHP 
DO 195 J=1»NHP 
195 C(I»J)=0.0 

DO 250 1=1»NHP 
DO 250 J=lrNSP 
DO 250 K=1»NHP 

250 C(I»J)=C(I»J)+B(I»K)*A(JrK) 

C 
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ooo ooo o o 


COMPUTE WD 

RS=RI+ (RO-RI ) * ( IM-2) / < IM-1 ) 

DO 500 1=1 * 1M 
DO 500 U=1 » KM 
K=IM*(J-1)+I 

R=RI+ (RO-RI ) * ( 1-1 ) / ( IM-1 ) 

TH=0 • 5235986* ( J“1 ) 

ARG=3.1415927*R/RS 
500 WD(K)=SIN(ARG>**2 

SUBTRACT OUT SPHERICAL PART FROM WD 

WB=0 • 0 

DO 300 IS=1#NSP 
JS=ISP(IS) 

I=ISP ( IS) /IM 
I=ISP(IS)-I*IM 

300 WB=WB+WD ( US ) ♦BE ( I ) 

DO 301 1=1 r IM 

301 Y(I)=BE(1)*WB/B2 
DO 302 IS=1 » NSP 
JS=ISP(IS) 

I=ISP(IS)/IM 

I=ISP(IS)-I*IM 

302 WD ( JS) =WD ( JS) -Y ( I ) 

IKM=IM*KM 

WRITE (6» 264) 

264 FORMAT (12H DISTURBANCE) 

WRITE (6# 265) (WD(K) »K=1» IKM) 

265 FORMAT (15E8» 3) 

DO 550 K=1 »NSP 
IS=ISP(K) 

550 WD (K) =WD (IS) 

COMPUTE PERFORMANCE INDEX 
AU1=0 . 0 

DO 600 1=1 » NSP 
600 AJ1=AJ1+WD ( I ) **2 
DO 700 1=1 » NSP 
DIJ=0.0 

DO 700 J=1»NSP 
IF(I.EQ.J) DI J=1 • 0 
DO 701 K=1 » NHP 
701 DI J=DI J-A (I»K)*C(K»J) 

700 WW(I)=WW(I)+DIU*WD(J) 

AJ2=0 • 0 

DO 750 1=1# NSP 
750 A J2=A J2+WD ( I ) *WW ( I ) 

WRITE (6# 1000 ) AJ1»AJ2 

1000 FORMAT (//46H THE PERFORMANCE INDEX BEFORE COMPENSATION IS »E10.5# 
1//45H THE PERFORMANCE INDEX AFTER COMPENSATION IS #E10.5) 

BJ1=SQRT ( AJ1/NSP) 

B J2=SQR T (AJ2/NSP) 

WRITE (6# 1050 ) BJ1»BU2 

1050 format ( 2 qh rms error before = #E2o.8/ 
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